Characteristics of Spatial–Temporal Variations in Coupling Coordination between Industrial Water Use and Industrial Green Development Systems in China

: The coupling coordination between industrial water use (IWU) and industrial green development (IGD) systems is crucial for achieving sustainable development goals. This paper measures the coupling coordination degree between IWU and IGD systems, and the spatial–temporal evolution characteristics of the degree are discussed. Here, the IWU system is assessed by adopting a dynamic slacks-based measure model, and the IGD system’s performance is evaluated using an entropy-weighted TOPSIS model. The results showed that: (1) The provincial IWU efﬁciency showed a rising trend from 2009 to 2018 in general, while its spatial dimension showed a distribution of high in the eastern region and low in the western region. (2) The IGD performance in the eastern region was better than that in the central and western regions; the gap in industrial innovation and industrial resources was the main factor. (3) The spatial structure of the coupling coordination degree between IWU and IGD systems was gradually stable. It also had a strong spatial dependence and its evolution volatility has been enhanced. (4) The coupling coordination was improving but exhibited a dynamic local spatial dependence and volatile process, and its spatial agglomeration had a relatively higher path dependence and locked spatial features


Introduction
Since 1978, rapid industrial progress has contributed irreplaceably to China's economic development [1,2]. However, due to the blind pursuit of economic development speed, rough industrial production leads to resource problems, and environmental management becomes the bottleneck restricting economic development. China's industry has yet to shake off the development pattern of high consumption and emissions [3], particularly evident in the industrial water use (IWU) system. In 2020, China's IWU reached 121.80 million tons, and the corresponding industrial wastewater discharge reached 2.07 million tons [4]. Among the direct pollution sources, the wastewater discharge from industrial pollution sources is second only to that of comprehensive sewage outlets [5], which is far behind the wastewater discharge standard of developed countries.
Since 2016, China has vigorously advocated the industrial green development (IGD) to promote high-quality industrial sustainable development in the 13th and 14th Five-Year Plans while protecting the environment and conserving resources. Made in China 2025 and Industrial Green Development Plan (2016-2020) pointed out that China's industrial development has entered the green transformation period. Besides, as a typical waterdeficient country, improving the IWU efficiency is particularly important in promoting IGD, which needs a series of policy guidelines and legal guarantees as the basis. The formulation of these policies and laws must be based on the scientific evaluation of IWU and IGD systems. Such assessment can also provide policymakers with quantitative information to help them conduct policy analysis effectively [6]. the temporal dimension for industrial adjustment and technological upgrading. However, traditional ESDA and mathematical statistics are difficult to represent either the characteristics spatial correlation or temporal dimension [39]. Exploratory spatial-temporal data analysis (ESTDA) can integrate temporal dimension analysis into spatial analysis more effectively and has been widely used in the analysis of land use, resource utilization, pollution control and other fields. Dong et al. [10] adopted ESTDA to reveal the spatial-temporal dynamics of the multidimensional poverty index in China from 2007 to 2017. By applying ESTDA, Zhang et al. [40] analyzed China's carbon footprint's spatial pattern and dynamic spatial-temporal interaction from 2001 to 2013. Liu et al. [41] used ESTDA and the spatial econometric model to analyze the spatial-temporal evolution characteristics and influencing factors of PM2.5 in the Yangtze River Economic Belt. Li et al. [42] analyzed the spatial-temporal dynamic evolution of GDP per capita for 2303 counties in China from 1998 to 2015 based on the ESTDA framework.
Based on the above analysis, we could find that existing research focuses more on the interaction between the IGD system and technological innovation or energy system, but few study the relationship between IWU and IGD systems. Besides, the method involved is not able to evaluate the dynamic evolution between systems effectively. In the coupling framework between resource and economic systems, this research focuses on the coupling coordination relationship between IWU and IGD systems and discusses its temporal and spatial variations to provide policy suggestions for achieving sustainable development goals. The objectives of this study are as follows: (1) Construct an index system to evaluate the performance of IWU and IGD systems. (2) Propose a new coupling coordination framework to evaluate the coupling coordination degree between IWU and IGD systems.
(3) Reveal the characteristics of spatial-temporal variations in coupling coordination degree by adopting ESTDA.
The remaining structure of the article is organized as follows. Section 2 introduces the construction of an index system, methods, and data acquisition. Section 3 presents the evaluation results and discussion of this research. In Section 4, the conclusions, limitations, and future research directions are provided.

Coupling Coordination Indicator and Methods
This paper constructs a coupling coordination framework between IWU and IGD systems, as shown in Figure 1. Firstly, A dynamic SBM model considering undesirable output referred by Tone [22] to evaluate the IWU efficiency from 2009 to 2018. Referring to existing research [19,43,44], we chose industrial employment, and water use as input indicators, industrial fixed assets as a carry-over indicator, industrial value added as desirable output indicators, and industrial wastewater discharge as undesirable output indicators.   The traditional TOPSIS model is difficult to eliminate the influence of the order of magnitude of index data on index weight. In order to eliminate this method bias, the entropy weight method is combined with the TOPSIS model. Based on the entropyweighted TOPSIS model referred by existing research [45,46], this paper puts forward a new index system to measure the performance of the IGD system, which involves industrial economy, industrial environment, industrial resources, and industrial innovation by integrating existing research [31][32][33] and Industrial Green Development Plan 2016-2020 issued by the Chinese government, as shown in Table 1. A coupling degree is a quantitative method for measuring the interaction in two or more systems. The coupling degree model of IWU and IGD systems is constructed by referring to relevant research [47], as shown in Equation (1).
where C is the coupling degree between IWU and IGD systems, 0 < C < 1, the greater the value of C is, the better the coupling is. f (x) and g(y) are the IWU efficiency and the performance of IGD, respectively. k is the adjustment coefficient, 2 ≤ k ≤ 5, which is determined by the number of subsystems in the coupling degree model, and the value of k is 2 in this research. The coupling coordination model is used to further discuss the coupling characteristics between IWU and IGD systems, and analyze the synergistic effect and overall efficacy of these two systems, as shown in Equation (2).
D is the coupling coordination degree of IWU and IGD systems, T is the comprehensive coordination index, α and β are the undetermined coefficients, and α + β = 1. Referring to previous research [48], this paper considers that IGD system is slightly more important than IWU system, α = 0.4, β = 0.6. The coupling coordination degree can be divided into five stages, as shown in Table 2.

Exploratory Spatial-Temporal Data Analysis
ESTDA integrates temporal and spatial features of geographic data to realize spatialtemporal interaction analysis. It can effectively reveal spatial distribution pattern changes of research objects [49,50]. ESTDA consists of LISA time path analysis, LISA spatialtemporal transition, spatial Markov chain, and other main research methods. LISA time path analysis and LISA spatial-temporal transition are adopted to reveal the characteristics of spatial-temporal variations in the coupling coordination degree.

LISA Time Path Analysis
LISA time path analysis presents the spatial-temporal cooperative changes and dynamic characteristics of geographic variables in the region by measuring the stability level of LISA coordinates composed of attribute variables and spatial lag in Moran's scatter plot [51]. The spatial-temporal interaction and dynamic characteristics of the coupling coordination degree between IWU and IGD systems are explained by the pairwise movement of the coupling coordination degree and its spatial lag value over time. The geometric characteristics of the LISA time path mainly include relative length (Γ i ) and tortuosity (f i ). Their respective expressions are shown in Equations (3) and (4).
L i,t is the LISA coordination of unit i in year t, and d(L i,t , L i,t+1 ) is the distance coordination of unit i from year t to year t + 1. If the coordination movement length of unit i is greater than the overall average level within the research time range, then Γ i is greater than 1. The larger Γ i is, the more significant the dynamic change of coupling coordination degree is, and vice versa. The larger f i is, the more curved LISA time path is. That is, unit i is more affected by neighborhood space (spillover/polarization) and its spatial evolution is more unstable.

LISA Spatial-Temporal Transition
LISA time path describes the geometric characteristics of the time trajectory of each space unit moving on Moran's I scatter plot, while the spatial-temporal transition further reveals the spatial relationship between the neighborhood in different local areas. Rey [50] embedded the attributes such as distance, direction, and condensation of each spatial unit . HH is the high-high state, indicating the positive synergistic movement/development of unit i and its adjacent units. LL is the low-low state, indicating the negative synergistic movement/development of unit i and its adjacent units. HL and LH are high-low and low-high states, respectively, representing the movement/development directions of unit i and its adjacent units are opposite.
There are four types of LISA spatial-temporal transitions: I, II, III and IV. Type I transition includes HH→HH, HL→HL, LH→LH, LL→LL, indicating that unit i and its adjacent units have no transition during the research period. Type II transition includes HH→LH, HL→LL, LH→HH, LL→HL, meaning that unit i has its transition, but its adjacent units do not have a transition. The type III transition includes HL→HH, HH→HL, LL→LH, LH→LL, indicating that unit i itself has no change, but its adjacent units have a transition. The type IV transition includes LL→HH, LH→HL, HL→LH, HH→LL, indicating that both unit i and its adjacent units have transitions. Further, type IV transition can be divided into types IV 1 (LL→HH, HH→LL) and IV 2 . The types IV 1 (LH→HL, HL→LH) indicate that unit i and its adjacent units have the same direction of transitions, while type IV 2 does the opposite. The details are shown in Table 3. Spatial-temporal flow (SF) and spatial-temporal condensation (SC) were used to characterize the path dependence of the spatial pattern of the coupling coordination degree between IWU and IGD systems, as shown in Equations (5) and (6).
In Equations (5) and (6), N j is the number of units in type j; M is the number of samples.

Research Area and Data Acquisition
The comparison object selected here is China's provincial administrative regions. By evaluating the relationship between IWU and IGD systems in each province, problems in industrial green transformation could be analyzed. In line with this idea, 30 provinces in China are selected as evaluation units, and the research period is from 2009 to 2018. Tibet, Hong Kong, Macau, and Taiwan are not included in the assessment due to the lack of data.
The total debt ratio of industrial enterprises above designated size, total profits of industrial enterprises above designated size, the full-time equivalent of industrial R&D employees, number of industrial R&D employees, number of industrial valid patents, and number of new industrial products are obtained from the database of the National Bureau of Statistics. Industrial tailpipe emissions, industrial wastewater discharge, industrial tailpipe purification, and industrial wastewater purification are brought from China Environmental Statistical Yearbook (2010-2020). Industrial employment, industrial water use consumption, industrial energy consumption, the industrial fixed assets stock, and the industrial value added are obtained from the China Industrial Statistical Yearbook (2010-2020) and China Statistical Yearbook (2010-2020). The remaining indicators are calculated from the corresponding raw data above.

Analysis of Industrial Water Use System
The IWU efficiency of each province during 2009-2018 is presented in Table 4. The IWU efficiency of most provinces was less than one, indicating that the IWU of these provinces needs to be optimized. On the whole, the IWU efficiency of each province presented an upward trend. From 2009 to 2018, the average IWU efficiency increased from 0.401 to 0.703, with an increase of about 75.19%. The improvement of IWU efficiency is closely related to the promulgation of water-related policies, including establishing a strict water resource management system, strictly controlling water for production, and effective wastewater treatment.  high water use consumption, adopted advanced water-saving technologies, and formulated strict water policies.
Ningxia, Gansu, Guizhou, Fujian, and Shaanxi were the five provinces with lower IWU efficiency. The mean values of their IWU efficiency were 0.409, 0.437, 0.441, 0.455, and 0.493, respectively. In these provinces, the problems, such as unbalanced economic development, low awareness of water conservation, and the widespread existence of industries with high water use consumption, lead to the waste of water resources in industrial production and low water use efficiency. However, there is no denying that the IWU in all the provinces mentioned above has improved since 2011 (the beginning of the 12th Five-Year Plan), which is due to the Chinese government's increasing emphasis on IGD, developing the recycling industry and optimizing the comprehensive utilization of water resources. Table 5 presents the index weights of China's IGD system based on the entropyweighted TOPSIS model. The four primary indexes indicators are ranked as follows: Industrial economy (27.23%) > Industrial innovation (26.51%) > Industrial resources (24.98%) > Industrial environment (21.28%). Economic development is still the primary goal of IGD, which is a necessary foundation to support the implementation of China's 14th Five-Year Plan. Industrial innovation is also a prerequisite required to improve industrial production efficiency, reduce resource consumption and reduce environmental pollution, and it plays the leading role in the IGD systems. The Chinese government has been increasing investment in industrial innovation in recent years. Industrial resource consumption and industrial pollution are important indexes to measure the sustainability and environmental friendliness of the IGD system. Table 5. The index weights of the IGD system.

Primary Indexes Weight Secondary Indexes Weight
Industrial economy 27.23% Labor productivity. X 1 7.30% Value-added growth rate. X 2 8.32% The total debt ratio of industrial enterprises above designated size. X 3 4.88% Total profits of industrial enterprises above designated size. X 4 6.73% Industrial environment 21.28% Industrial tailpipe emission. X 5 3.83% Industrial wastewater discharge. X 6 5.32% Industrial tailpipe purification. X 7 5.11% Industrial wastewater purification. X 8 2.32% The proportion of operating expenses of wastewater purification facilities in value added.

2.50%
The proportion of operating expenses of tailpipe purification facilities in value added.

2.20%
Industrial resources 24.98% Industrial energy consumption. X 11 7.33% Industrial water use consumption. X 12 7.21% Industrial energy intensity. X 13 5.55% Industrial water use intensity. X 14 4.89% Industrial innovation 26.51% Full-time equivalent (FTE) of industrial R&D employees. X 15 7.14% Number of industrial R&D employees. X 16 8.70% Number of industrial valid patents. X 17 5.30% Number of new industrial products. X 18 5.37% The top five secondary indexes are as follows: Number of industrial R&D employees (8.70%) > value-added growth rate (8.32%) > industrial energy consumption (7.33%) > labor productivity (7.30%) > industrial water use consumption (7.21%). Number of industrial R&D employees is the industrial innovation index. China's industrial structure is transforming; energy-intensive and inefficient heavy industries still occupy a large proportion of the industry, which is a significant cause of high energy and resource consumption and heavy pollution. Increasing investment in industrial innovation can not only support the adjustment of industrial structure but also drive the IGD [51]. Value-added growth rate and Labor productivity are the important indexes to measure the performance of the industrial economy. As the supporting industry of the national economy, industry bears the critical responsibility of promoting social and economic development. The other two indexes are industrial resources indexes. Resource consumption and environmental pollution is the bottleneck restricting the IGD. Since the 12th Five-Year Plan period, the Chinese government has made great efforts to introduce advanced industrial production technology and increase investment in environmental pollution control to ensure the green transformation of industrial enterprises. Table 6 shows the performance of the IGD system of 30 provinces in China. From 2009 to 2018, IGD has been improved in most provinces. Zhejiang, Jiangsu, Shanghai, Tianjin, Beijing, and Guangdong were the six provinces with outstanding performance in the IGD system, and the mean values of their performance of the IGD system exceeded 0.4. They are all the most advanced industrial provinces in China and are all located in the eastern region. Compared with other provinces, they have a complete and efficient industrial structure and are the earliest provinces to promote industrial adjustment. Based on a strong economy, they have vigorously developed high-tech industries, increased investment in environmental pollution control, and eliminated industries with high resource consumption and pollution. Excellent transportation and geographical conditions also ensure their access to advanced production technology, foreign investment, and policy preferences. The provinces in the central region, such as Anhui, Shanxi, and Jilin, were in the second group. Although their performance was not as good as that of the provinces in the eastern region, it should not be ignored that their IGD has maintained an annual growth rate of more than 4%. In recent years, the provinces in the central region have increased cooperation with provinces in the eastern region to introduce advanced energy-saving and water-saving technologies. Since the implementation of the Plan on the Rise of the Central Region [52], the Chinese government has increased the policy preference for the provinces in the central region, especially allocating funds to recruit R&D employees and attract investment for these provinces, to build them into modern equipment manufacturing and high-tech industrial bases. At the same time, the government strictly controlled environmental pollution to ensure the sustainable development of the industrial economy.
The provinces in the western region were at the bottom in terms of the performance of the IGD system during 2009-2018. As the only municipality directly under the central government in the western region, the industrial development of Chongqing was very rapid. Except for Chongqing, the IGD degree of other provinces in the western region increased slowly, and even Guangxi, Ningxia and Qinghai showed a downward trend. Compared with other provinces in the western region, Chongqing is more attractive to foreign investment and has easier access to high technology, supporting its IGD. Other provinces in the western region, in order to ensure economic growth, undertake the backward industries eliminated from the central and eastern regions more, resulting in the intensification of resource consumption and environmental pollution.
As shown in Table 6 and Figure 2, there were only two provinces (Jiangsu and Zhejiang) whose performance of the IGD system was greater than 0.4 in 2009, while in 2018, there were six provinces (Jiangsu, Zhejiang, Shanghai, Beijing, Tianjin, Shandong and Chongqing) whose IGD degree was greater than 0.4. Figure 2 also presents the composition of the IGD degree in 2009 and 2018. Obviously, except for Shaanxi, the industrial economic index score ratio to the IGD degree has decreased in the other 29 provinces. The slowdown in economic growth is one reason for this phenomenon. More importantly, the Chinese government is paying more and more attention to water consumption and industrial environment control instead of focusing solely on the industrial economy.

Analysis of Coupling Coordination
The coupling coordination model is used to calculate the coupling coordination degree between IWU and IGD systems

Analysis of Coupling Coordination
The coupling coordination model is used to calculate the coupling coordination degree between IWU and IGD systems during 2009-2018, and the results are shown in During the research period, the coupling degree of the provinces in the eastern region remained high, and the fluctuation range was small. Their coupling degree was mainly in the high coupling stage and run-in stage. The coupling degree of provinces in the central region decreased and then increased. The coupling development between IWU and IGD systems in these provinces experienced a process from deterioration to improvement. The coupling degree of provinces in the western region was generally low. Except for Chongqing, the coupling degree of other provinces was mainly in low couple stage and irrelevant stage. As shown in Figure 4a     In these provinces, the construction of efficient IWU system and wellconsidered IGD system ensure that their coupling coordination degree has been maintained at a high value. On the contrary, the coupling coordination degree of the other 7 provinces (Guangxi, Hainan, Liaoning, Gansu, Qinghai, Ningxia, and Xinjiang) all decreased during the research period. Among these seven provinces, the IWU efficiency of Liaoning and Hainan has been effectively improved. Nonetheless, their IGD degree was always low, which affected the coupling and coordinated development between IWU and IGD systems. Guangxi, Gansu, Ningxia, Xinjiang, and Qinghai performed poorly in the IWU and IGD systems, ranking at the bottom of the 30 provinces. These provinces scored low in industrial resources, industrial environment, and innovation, leading to the decline of their coupling coordination degree between the IWU and IGD systems.
Obviously, as shown in Figure 4, the research period can be divided into two stages: 2009 to 2013 ( Figure 4a) and 2014 to 2018 (Figure 4b). During 2009-2013, 30 provinces have developed just four kinds of coupling coordination stages. The numbers of provinces in the high dysregulation stage were 0, 0, 5, 5, and 3, respectively. The numbers of provinces in the basic dysregulation stage were 11, 13, 12, 11, and 14, respectively. The numbers of provinces in the run-in stage were 15, 13, 7, 6, and 6, respectively. The numbers of provinces in the basic coordination were 4, 4, 6, 8, and 7, respectively. During this period, the coupling and coordinated development between IWU and IGD systems deteriorated. There were no provinces in the high coordination stage, the number of provinces in the high dysregulation stage and basic dysregulation stage showed an upward trend. Provinces in the run-in stage were gradually changed to basic dysregulation stage and high dysregulation stage. During 2014-2018, 30 provinces developed five kinds of coupling coordination stages. The numbers of provinces in the high dysregulation stage were 0, 3, 1, 3, and 0, respectively. The number of provinces in the basic dysregulation stage was 9, 11, 13, 10, and 8, respectively. The number of provinces in the run-in stage was 11, 7, 7, 8, and 10, respectively. The provinces in the basic coordination stage were 7, 6, 5, 8, and 9, respectively. The provinces in the high coordination stage were 7, 6, 5, 8, and 9, respectively. Compared with the previous period, the provincial coupling coordinated development has been improved. Number of provinces in the high dysregulation stage and basic dysregulation stage showed a downward trend, the number of provinces in the high coordination stage increased gradually. Provinces in the run-in stage and basic dysregulation stage gradually changed to the basic coordination stage. In these provinces, the construction of efficient IWU system and well-considered IGD system ensure that their coupling coordination degree has been maintained at a high value. On the contrary, the coupling coordination degree of the other 7 provinces (Guangxi, Hainan, Liaoning, Gansu, Qinghai, Ningxia, and Xinjiang) all decreased during the research period. Among these seven provinces, the IWU efficiency of Liaoning and Hainan has been effectively improved. Nonetheless, their IGD degree was always low, which affected the coupling and coordinated development between IWU and IGD systems. Guangxi, Gansu, Ningxia, Xinjiang, and Qinghai performed poorly in the IWU and IGD systems, ranking at the bottom of the 30 provinces. These provinces scored low in industrial resources, industrial environment, and innovation, leading to the decline of their coupling coordination degree between the IWU and IGD systems.
Obviously, as shown in Figure 4, the research period can be divided into two stages: 2009 to 2013 ( Figure 4a) and 2014 to 2018 (Figure 4b). During 2009-2013, 30 provinces have developed just four kinds of coupling coordination stages. The numbers of provinces in the high dysregulation stage were 0, 0, 5, 5, and 3, respectively. The numbers of provinces in the basic dysregulation stage were 11, 13, 12, 11, and 14, respectively. The numbers of provinces in the run-in stage were 15, 13, 7, 6, and 6, respectively. The numbers of provinces in the basic coordination were 4, 4, 6, 8, and 7, respectively. During this period, the coupling and coordinated development between IWU and IGD systems deteriorated. There were no provinces in the high coordination stage, the number of provinces in the high dysregulation stage and basic dysregulation stage showed an upward trend. Provinces in the run-in stage were gradually changed to basic dysregulation stage and high dysregulation stage. During 2014-2018, 30 provinces developed five kinds of coupling coordination stages. The numbers of provinces in the high dysregulation stage were 0, 3, 1, 3, and 0, respectively. The number of provinces in the basic dysregulation stage was 9, 11, 13, 10, and 8, respectively. The number of provinces in the run-in stage was 11, 7, 7, 8, and 10, respectively. The provinces in the basic coordination stage were 7, 6, 5, 8, and 9, respectively. The provinces in the high coordination stage were 7, 6, 5, 8, and 9, respectively. Compared with the previous period, the provincial coupling coordinated development has been improved. Number of provinces in the high dysregulation stage and basic dysregulation stage showed a downward trend, the number of provinces in the high coordination stage increased gradually. Provinces in the run-in stage and basic dysregulation stage gradually changed to the basic coordination stage.

LISA Time Path Analysis
Based on the evolution characteristics of the coupling coordination degree shown in Section 3.3 (Figures 3 and 4), the research period was divided into two development stages: 2009-2013 and 2014-2018. The relative length, tortuosity, and movement direction of coupling coordination degree between IWU and IGD systems in China were calculated based on the LISA time path, as shown in Figure 5.
Based on the natural break point method in ArcGIS and manual classifications, the relative lengths of the LISA time path were classified into a low relative length (Γ i < 0.33), medium relative length (0.33 < Γ i < 0.9), higher relative length (0.9 < Γ i < 1.3), and high relative length (Γ i > 1.3) (Figure 5a,b). During 2009-2013, there were 11 provinces with a relative length greater than one. There were nine provinces with higher relative length and six provinces with high relative length, these provinces included Xinjiang, Yunnan, Jilin, Shanxi, Henan, and Anhui, indicating that the local spatial structure was highly dynamic. In contrast, there were 14 provinces with medium relative length, and only Shanghai was in low relative length. These provinces, such as Shanghai, Jiangsu, Shandong, Fujian, Anhui, Beijing, Tianjin, Chongqing, and Sichuan, are mostly located in the eastern and central region, with stable local spatial structures. During 2014-2018, there were 14 provinces with a relative length greater than 1. Although the proportion of provinces with relative length greater than the national average increased, only Henan was in high relative length. Moreover, the number of provinces with low relative length increased to seven. Generally, during 2009-2018, the spatial structure of the coupling coordination degree between IWU and IGD systems was gradually stable and the spatial structure of the eastern region is more robust than that of the central and western regions.
Similarly, the tortuosity of the LISA time path was classified into a low tortuosity (0.45 < f i < 3.24), medium tortuosity (3.24 < f i < 8.22), higher tortuosity (8.22 < f i < 19.58), and high tortuosity (19.58 < f i < 67.32) (Figure 5c,d). During 2009-2013, In all provinces except Shanghai, the tortuosity was more significant than 1, indicating that most provinces were affected by the spatial dependence of local structures to varying degrees. Gansu and Henan were in high tortuosity, and Jiangsu and Qinghai were in higher tortuosity. It reflected the solid dynamic characteristics of their spatial transition process. A total of 18 provinces were in low tortuosity, accounting for 60%, and these provinces presented a continuous patchy distribution. During 2014-2018, the tortuosity of all provinces was more significant than 1, indicating that all provinces were affected by the spatial dependence of local structures. Jiangsu and Zhejiang were in high tortuosity, Guangdong, Qinghai, Gansu, Inner Mongolia, Hebei, Shandong, and Henan were in higher tortuosity, thus reflecting the intense volatility of the spatial dependence of coupling coordination degree in these provinces and a dynamic change process among these provinces and their neighborhoods. The provinces with low tortuosity decreased to only seven, accounting for 23.33%. In contrast, the proportion of provinces with high and higher tortuosity increased to 30%, indicating that the coupling coordination degree between IWU and IGD systems had a strong spatial dependence. Its evolution volatility has been enhanced during the research period.
The LISA movement direction was used to characterize the spatial integration of the spatial variation in the coupling coordination between IWU and IGD systems in China. Four types were delineated by referring to the average level. As shown in Figure 5e,f, during 2009-2013, Jiangsu, Zhejiang, Shanghai, Beijing, Tianjin, Shandong, Fujian, Henan, and Hubei showed a synergistic high growth trend, while Chongqing, Guizhou, Guangdong, and Shanxi showed a synergistic negative growth trend. Provinces with high synergistic growth are mainly located in the eastern region, while those with low synergistic growth are primarily located in the western region. During 2014-2018, the number of provinces with high synergistic growth decreased from eight to seven, and the number of provinces with low synergistic growth increased from three to five, indicating that the spatial integration of the spatial evolution in coupling coordination was slightly enhanced in this period compared to the previous period. The spatial variation of the coupling coordination degree between IWU and IGD systems in China performed differently in the east, central, and western regions. This situation remains consistent in the short term.

LISA Time Path Analysis
Based on the evolution characteristics of the coupling coordination degree shown in Section 3.3 (Figures 3 and 4), the research period was divided into two development stages: 2009-2013 and 2014-2018. The relative length, tortuosity, and movement direction of coupling coordination degree between IWU and IGD systems in China were calculated based on the LISA time path, as shown in Figure 5. Based on the natural break point method in ArcGIS and manual classifications, the relative lengths of the LISA time path were classified into a low relative length ( Γ i < 0.33), medium relative length (0.33 < Γ i < 0.9), higher relative length (0.9 < Γ i < 1.3), and high relative length ( Γ i > 1.3) (Figure 5a,b). During 2009-2013, there were 11 provinces with a relative length greater than one. There were nine provinces with higher relative length and six provinces with high relative length, these provinces included Xinjiang, Yunnan,

LISA Spatial-Temporal Transition
This paper constructs the spatial-temporal transition matrix to explore the transfer characteristics of the local correlation types of coupling coordination degree between IWU and IGD systems. In Table 7, from 2009 to 2018, LL/LL had the highest transfer probability (0.333) with 91 samples. In terms of the transition type, the most common transition was Type I, which had a ratio of 82%, thus indicating that the spatial agglomeration of coupling coordination degree between IWU and IGD systems had a relatively higher path dependence and locked spatial features. The reason was that, in the process of rapid industrialization, the low IWU efficiency and low IGD degree caused by the resourceintensive industrial development mode could not be changed in the short term. The ratios of transition types II, III, and IV were 0.078, 0.067, and 0.033, respectively, thus reflecting that the internal factors of the unit had a more significant impact on their concentration of coupling coordination than the neighborhood overflow. The SF was only 0.144, and SC was 0.841 from 2009 to 2018. It indicates that the distribution of coupling coordination degree between IWU and IGD systems in China has strong transfer sluggishness, and the coupling coordination degree was relatively stable.

Conclusions
This paper used the dynamic SBM model and the entropy-weighted TOPSIS model to measure the IWU efficiency and the performance of the IGD system of 30 provinces in China, respectively. Further, a new coupling coordination framework measures the interaction between IWU and IGD systems. The spatial-temporal pattern and path evolution characteristics of the coupling coordination degree are reasonably analyzed using the ESTDA framework. The main conclusions are as follows: (1) In terms of temporal dimension, the IWU efficiency of each province showed an overall rising trend. In contrast, in terms of spatial extent, it showed a distribution of high in the east region and low in the west region. The IGD degree of the provinces in the eastern region was better than that in the central and western regions; the main reason is the gap in industrial innovation and resources. (2) The coupling coordination between IWU and IGD systems was getting better. Overall, the number of provinces in high coordination showed an increasing trend, and the number of provinces in high dysregulation showed a downward trend. Nevertheless, the coupling coordination degree showed apparent regional differences in the eastern, central, and western regions, and some poorly performing provinces still need to improve efficiently. (3) The spatial structure of the coupling coordination degree between IWU and IGD systems was gradually stable, and the spatial structure of the eastern region is more robust than that of the central and western regions. The coupling coordination degree between IWU and IGD systems had a strong spatial dependence, and its evolution volatility has been enhanced during the research period. (4) There was a dynamic local spatial dependence and volatile process in the coupling coordination degree between IWU and IGD systems. Its spatial agglomeration had a relatively higher path dependence and locked spatial features.
Based on the above results, some relevant policy suggestions are put forward. First, strengthen regional governance, and formulate "common but differentiated" industrial water resources utilization and IGD strategy. The eastern region should lead in efficiently utilizing IWU and IGD systems. The central and western regions should close down industries with high water use consumption and reduce excess capacity while actively protecting the environment and increasing investment in innovation. Second, establish the regional joint prevention and control system to supervise industrial water conservation and environmental protection. Third, vigorously develop green industries, adjust industrial structure, and advocate low water consumption industrial economy to achieve industrial sustainable development.

Limitations and Future Research Directions
By constructing the coupling coordination framework between IWU and IGD systems, this paper provides a reference for the coordinated development of ecology and economy, but there are still areas for improvement to be filled. Due to the limited data availability, we did not include research data after 2020 in this research. However, it is undeniable that after 2020, China's industrial green development has been rapid, which has a very representative research value. This also leads to a decrease in the representativeness of the research. In terms of pollutant discharge, the pollutant discharge considered in this paper mainly includes wastewater discharge and tailpipe emission. We did not consider solid waste discharges that might influence the research results. At the same time, different types of tailpipe emission and wastewater discharge have a massive difference in the impact of environmental pollution and management, which is not included in the analysis. The influence of government policy priorities and residents' environmental awareness is also essential. The data and content must be further supplemented and improved in future research. The factors that affect the coupling coordination between IWU and IGD systems are complex and affected by complex factors at different scales. More is needed to discuss the spatial-temporal distribution of the coupling coordination degree. To ensure the reliability and universality of the research, the influence mechanism of various factors on coupling coordination degree should be discussed emphatically in future research.