Regional Coordination and Security of Water–Energy–Food Symbiosis in Northeastern China

The interactions of water, energy, and food, which are essential resources for human survival, livelihoods, production, and development, constitute a water–energy–food (WEF) nexus. Applying symbiosis theory, the economic, social, and natural factors were considered at the same time in the WEF system, and we conducted a micro-level investigation focusing on the stability, coordination, and sustainability of the symbiotic units (water, energy, and food), and external environment of the WEF system in 36 prefecture-level cities across three northeastern provinces of China. Finally, we analyzed the synergistic safety and coupling coordination degree of the WEF system by the combination of stability, coordination, and sustainability, attending to the coordination relationship and influences of the external environment. The results indicated that the synergistic safety of the WEF system in three northeastern provinces need to equally pay attention to the stability, coordination, and sustainability of the WEF system, since their weights were 0.32, 0.36 and 0.32, respectively. During 2010–2016, the synergistic safety indexes of the WEF system ranged between 0.40 and 0.60, which was a state of boundary safety. In the current study, the coupling coordination degree of the WEF system fluctuated around a value of 0.6, maintaining a primary coordination level; while in the future of 2021–2026, it will decline to 0.57–0.60, dropping to a weak coordinated level. The conclusion could provide effective information for decision-makers to take suitable measures for the security development of a WEF system.


Introduction
The water-energy-food (WEF) nexus is a concept that describes the complex underlying mechanisms behind mutually constraining or reinforcing interactions of water, energy, and food in the current global context of resource scarcity. As a result of population growth and rapid economic and social development, demands for water, energy, and food are rising and are expected to increase by around 40%, 50%, and 35%, respectively, by 2030 [1]. Concurrent population expansion, resource shortages, and environmental degradation critically affect the WEF system, and meeting current demands for food, energy, and water for societal development poses a major challenge [2]; for example, the changing climate is expected to aggravate water and energy securities [3].
Water, energy, and food have complex interconnections. Water is required to produce energy, while energy is needed for water extraction, treatment, and distribution. The food sector requires water and energy to produce food products, while fertilizer and pesticides from farmland have a negative impact on water quality; however, biomass is a potential alternative energy source [4]. The fluctuations in any one subsystem will affect the others [5]. From an integrated multi-resource perspective, conducting research on the WEF nexus, strengthening inter-departmental coordination and cooperation, and safeguarding regional water, energy, and food security have become top concerns for the international community [6]. In January 2011, the World Economic Forum released its Global Risks ronmental impact controls to optimize regional water, energy, and food supply strategies. Bai et al. analyzed changes in pressures relating to the supply and demand of water, energy, and food in various regions of China [36]. However, previous studies mainly focused on WEF system linkages at national and provincial levels [37], linkage frameworks [38], measurement of safety risks [39], and optimized synergistic development [40]. Few studies focus on the spatial and temporal changes of WEF linkages at local and municipal levels. The indicators used to evaluate the security of the WEF system in previous researches mainly focused on three dimensions: water resources, energy, and food, however, with the external environment (social, economic, and natural) receiving less attention.
A comprehensive theoretical approach is necessary to apply and detect the reciprocal process of water, energy, and food elements in the WEF system. Recently, the symbiosis theory was introduced in the research of WEF and was helpful for better understanding the complex relationship between WEF and its external environment [41]. In the context of the WEF system, symbiotic relationships are fostered among system elements through coupling mechanisms that ensure the stable and harmonious development of the system [42]. Based on the concept of symbiosis [43,44], these three elements (water, energy, and food) can be considered as the symbiotic units in the development of the WEF system. These units interact directly at the symbiotic interface. At the same time, changes in the social, economic, and ecological environments of the WEF system affect these interactions. Based on the symbiosis theory, Zhi et al. [41] analyzed the WEF nexus in the whole of China and applied the Press-State-Impact-Response (PSIR) model to construct the assessment index system of fitness. The results showed the fitness level in northeast and eastern China were relatively high while the fitness number in the central and western China were relatively low.
The northeastern China consists of three large industrial and commercial grainproducing provinces; the synergistic safety level of WEF in this region is crucial to the stable development of the total system. In this study, based on the symbiosis theory, we considered the economic, social, and natural factors in the WEF system. By the combination of comprehensive evaluation method and a coupling model, we evaluated the levels of synergistic safety and coupling coordination of the WEF systems in 36 prefecture-level cities across the three provinces in northeastern China from 2010 to 2016 (see Supplementary file). This study has three objectives: (1) to construct a safety evaluation system for WEF systems and to evaluate the synergistic safety level of each prefecture-level municipality; (2) to predict coupling coordination degrees in the three northeastern provinces using the gray model (GM) during the period 2017-2026; and (3) to explore the main impact indicators of the WEF nexus in the three provinces and to hopefully provide scientific measures for guiding the sustainable development of the social, economic, and ecological environments within these provinces.

Study Area
The study area encompasses 36 prefecture-level cities across three northeastern Chinese provinces (38 • 43 -53 • 33 N, 118 • 53 -135 • 05 E): Heilongjiang, Liaoning, and Jilin Provinces (Figure 1), covering an area of approximately 787,300 km 2 and accounting for about 8.2% of China's total land area. The Songnen Plain, Sanjiang Plain, and Liaohe Plain are three of the principal regions of agricultural production. In 2018, the grain-producing region of the three northeastern provinces reached 13.33 million tons, accounting for about 20.3% of the total amount in China. In 2017, the total water consumption in Heilongjiang Province, Jilin Province, and Liaoning Province was 7.43 × 10 10 m 3 , 3.94 × 10 10 m 3 , and 1.86 × 10 10 m 3 , respectively. The agricultural water consumption accounted for 89.6%, 66.1%, and 53%, respectively, and the industrial water consumption accounted for 5.6%, 4.3%, and 13.6%, respectively. The main oil fields in the three northeast provinces are Daqing, Jilin and Liaohe. Daqing Oilfield is China's largest crude oil supply base. Liaohe Oilfield is the third largest, with oil and natural gas reserves accounting for 15% and 10% of China's oil, respectively. In 2019, the gross domestic product (GDP) of the three northeastern provinces was 502.46 billion renminbi (RMB), accounting for 5.07% of the total national.
Province, Jilin Province, and Liaoning Province was 7.43 × 10 10 m 3 , 3.94 × 10 10 m 3 , and 1.86 × 10 10 m 3 , respectively. The agricultural water consumption accounted for 89.6%, 66.1%, and 53%, respectively, and the industrial water consumption accounted for 5.6%, 4.3%, and 13.6%, respectively. The main oil fields in the three northeast provinces are Daqing, Jilin and Liaohe. Daqing Oilfield is China's largest crude oil supply base. Liaohe Oilfield is the third largest, with oil and natural gas reserves accounting for 15% and 10% of China's oil, respectively. In 2019, the gross domestic product (GDP) of the three northeastern provinces was 502.46 billion renminbi (RMB), accounting for 5.07% of the total national.

Data
Four types of data were used in this study, namely: food, energy, water resources, and social economic statistical data. The food data included total food production (rice, wheat, maize, sorghum, cereals, potatoes, and soybeans), sown area, and arable land area, which were obtained from statistical yearbooks. The energy data, including energy production, energy consumption (primary, secondary, and tertiary industries as well as residential energy consumption), and total power of agricultural machinery, were obtained from statistical yearbooks (https://www.epsnet.com.cn/index.html#/Home;https://data.cnki.net/Yearbook/Navi?type=type&code=A). Given the unavailability of data on energy production and consumption in some prefecture-level cities, we based our calculations on data on the outputs of provincial energy production and energy consumption as a proportion of the GDP of various cities, that is, = • , with X denoting each prefecture-level city and T denoting the three northeastern provinces [45]. The water used data including water supply, total water resources, agricultural water consumption, and total water consumption were obtained from water resources bulletins (http://www.mwr.gov.cn/sj/tjgb/szygb/;

Data
Four types of data were used in this study, namely: food, energy, water resources, and social economic statistical data. The food data included total food production (rice, wheat, maize, sorghum, cereals, potatoes, and soybeans), sown area, and arable land area, which were obtained from statistical yearbooks. The energy data, including energy production, energy consumption (primary, secondary, and tertiary industries as well as residential energy consumption), and total power of agricultural machinery, were obtained from statistical yearbooks (https://www.epsnet.com.cn/index.html#/Home; https://data. cnki.net/Yearbook/Navi?type=type&code=A). Given the unavailability of data on energy production and consumption in some prefecture-level cities, we based our calculations on data on the outputs of provincial energy production and energy consumption as a proportion of the GDP of various cities, that is, X Energy Production = X GDP T GDP ·T Energy Production , with X denoting each prefecture-level city and T denoting the three northeastern provinces [45]. The water used data including water supply, total water resources, agricultural water consumption, and total water consumption were obtained from water resources bulletins (http://www.mwr.gov.cn/sj/tjgb/szygb/; http://slt.ln.gov.cn/jbgb/szygb/; http: //slt.jl.gov.cn/zwgk/szygb/) published for each of the three northeastern provinces and the prefecture-level cities. The social economic statistical data including total population, GDP, and fertilizer application were obtained from statistical yearbooks (https: //data.cnki.net/Yearbook/Single/N2019040178).
As the study area encompassed 36 prefecture-level cities and their data have only been updated up to 2016 in statistical yearbooks, on the basis of ensuring the completeness of the data, we chose the period 2010-2016. The intrinsic parameters are the factors reflecting the intrinsic properties of symbiotic units. There should be at least one set of intrinsic parameters that are compatible across the symbiotic units. The principal intrinsic parameter of the water, energy, and food subsystems within the WEF system of the three northeastern provinces is the productive capacity of the food subsystem. The WEF system are heterogeneous because of different symbiotic units. Increased food production has driven the formation of the symbiotic interfaces of the WEF system in the three northeastern provinces, whereas declining water resources and increased energy consumption have been opposing forces. As food production expands, water and energy resources need to be conserved and maintained according to their carrying capacities. A stable relationship among the subsystems is maintained and principal intrinsic parameter of the symbiotic units remains compatible through existing linkages among food cultivation, development, and processing of water and energy resources. When a symbiotic relationship is absent, that is, when food production is characterized by predatory exploitation of water and energy resources, depletion of water resources and degradation of the environment are inevitable outcomes, ultimately leading to a situation in which food production is no longer possible. The utilization of water and energy resources should therefore be adjusted in the course of food production to foster symbiotic relations among the subsystems of the WEF system. Water, energy, and food subsystems interact under the socio-economic-natural system to achieve the overall development of food, that is, to generate symbiotic energy [41,46,47] (Figure 2).

Construction of the Evaluation Index System
We developed an index system for evaluating the synergistic security of the WEF system in the three northeastern Chinese provinces according to the principles of systematicity and comprehensiveness, scientificity and practicality, regionality and comparability, and operability and quantification [48]. We constructed three indexes for evaluating symbiosis: a stability index, a coordination index, and a sustainability index. The stability

Construction of the Evaluation Index System
We developed an index system for evaluating the synergistic security of the WEF system in the three northeastern Chinese provinces according to the principles of systematicity and comprehensiveness, scientificity and practicality, regionality and comparability, and operability and quantification [48]. We constructed three indexes for evaluating symbiosis: a stability index, a coordination index, and a sustainability index. The stability index measured the quantity, quality, structure, function, and carrying capacity of the water, energy, and food subsystems. The main indicators were the rates of development and utilization of water resources and water consumption per 10,000 RMB of GDP, energy self-sufficiency and consumption per 10,000 RMB of GDP, and the rate of utilization of cultivated land and food production per capita. Because the three northeastern provinces are major food-producing provinces in China, the coordination index measured resource allocation and utilization efficiency in the process of transforming water-food and energyfood dependence. The main indicators were the proportion of agricultural water to the total land, the amount of water use per unit of food production, the power consumption of agricultural machinery per unit of cultivated land, and the proportion of energy consumption within the primary industrial sector. The sustainability index measured the degrees of adaptation and interaction among the WEF system and the external environment using the indicators of GDP per capita and chemical fertilizer consumption per unit of sown land [49][50][51][52][53][54] (Table 1). Based on the index system, we applied a comprehensive evaluation method and a coupling model which were introduced in Sections 2.3.2 and 2.3.3 to establish the synergistic safety index and coupling coordination index, then to evaluate the levels of synergistic safety and coupling coordination of the WEF system, and to predict the coupling coordination in the period of 2017-2026 ( Figure 3).

Entropy Method
The entropy assignment method is an objective method applied to determine the weights of indicators according to the amount of information that they provide, thus effectively avoiding deviations attributed to subjective factors [55]. Nowadays, this method has been widely applied to various research fields to determine factor weights [56][57][58]. The first step in the calculation entailed processing the original dimensionless data to eliminate the complexity associated with variations in the different dimensions of the data and to achieve comparability and measurability of the indicators. The following weighting process was applied: Positive indicators: Negative indicators:

Entropy Method
The entropy assignment method is an objective method applied to determine the weights of indicators according to the amount of information that they provide, thus effectively avoiding deviations attributed to subjective factors [55]. Nowadays, this method has been widely applied to various research fields to determine factor weights [56][57][58]. The first step in the calculation entailed processing the original dimensionless data to eliminate the complexity associated with variations in the different dimensions of the data and to achieve comparability and measurability of the indicators. The following weighting process was applied: Positive indicators: Negative indicators: where X ij is the standardized indicator value. X ij denotes the value of the jth indicator of the ith sample, and its weight (P ij ) was calculated as follows: The entropy value for indicator j was calculated using the following equation: where k > 0, e j > 0, and k is related to the sample size m. In general, if k = 1/ln m , then 0 ≤ e ≤ 1. The coefficient of variation for indicator j was calculated as follows: Greater variation in the indicator value X ij corresponded to a smaller entropy value. Moreover, a higher value of g i indicated that it had a more significant role in the overall evaluation. The weights of the indicators were conducted as follows:

Comprehensive Evaluation Method
Using the constructed index system for evaluating the WEF nexus and the determined weights, the following equations were obtained in a comprehensive evaluation [59] of the WEF system in the three northeastern provinces: where S(x), C(y), and E(z) are the stability, coordination, and sustainability indices; a, b, and c are the respective weights of each indicator, x , y , and z are standardized data; and i, j, and k are the numbers of indicators selected within each system. Based on the stability, coordination, and sustainability indices, we calculated the level of synergistic safety of the WEF system as follows: where T is the synergistic safety level of the WEF system; and α, β, and γ are the weights of influence of each subsystem on social development. In this study, the synergistic safety score of the WEF system ranged between 0 and 1, with a step size of 0.2. We distinguished five safety levels: extremely unsafe, unsafe, boundary safe, safer, and very safe.

Construction of the Coupling Model
The degree of coupling indicates the extent to which systems or elements interact with each other, while the degree of coordination indicates whether the system is functioning well. We constructed the following coupling model on the basis of a review of relevant studies [60]: Incorporating the three subsystems (n = 3), namely the stability, coordination, and sustainability subsystems, we calculated the coupling degree using the following equation: where C denotes the coupling degree and C ∈ [0, 1]. Referring to the relevant literature and to the principles of the coupled development of WEF system in the study area [60], we divided the coupling degree into four levels: low-level coupling (C ∈ [0, 0.3]), coupling (C ∈ (0.3, 0.5]), run-in (C ∈ (0.5, 0.8]), and high-level coupling indicating the commencement of the stage of optimized coupling (C ∈ (0.8, 1.0]). Having calculated the coupling degree (C ) and the synergistic safety level (T), we calculated the coupling coordination degree (D) as follows: The uniform distribution function method was used to determine the type of coupling coordination and the classification criteria [60], as shown in Table 2. The GM (1, 1) model consists of a single variable's first-order component equation and can be used to resolve the issue of uncertain prediction problems with at least four numbers. The use of differential equations can be applied to fully explore the essence of information to achieve predictions that are highly accurate. Irregular raw data can be organized into regular sequences [61].
If the time series X (0) = X (0) (1), X (0) (2), . . . , X (0) (n) has n observations, and a new series X (1) = X (1) (1), X (1) (2), . . . , X (1) (n) is generated through a process of information accumulation, then the corresponding differential equation for the GM (1, 1) model can be calculated using the following equations: where X (1) is the new sequence generated through the accumulation of n sequence values, t is the nth sequence value, a is the development coefficient, and b is the amount of the grey effect. The parameter vectors can then be solved using the least squares method: where B denotes the data matrix; Y is the data vector, and the predictive model is expressed as:

Model Testing
Residual, correlation, and post-hoc tests are commonly performed for grey prediction. However, in this study we performed the posterior error test, which entails the following calculations: (1) The standard deviation of the original sequence (S 1 ): (2) The standard deviation of the absolute error series (S 2 ): (3) The variance ratio (C): (4) The small error probability (P): Order: Then P = P{e i < S 0 }. Table 3 the criteria used to categorize the levels of model validation [61].

Impact Indicators
The weights of the evaluation indicators reflect the degree of influence that each indicator has on the security of the WEF system. As shown in Figure 4, during the period 2010-2016, the indicators that had the strongest influence on the security of the WEF system in the three northeastern provinces, in descending order of importance, were: GDP per capita, agricultural water use as a proportion of total water use, energy self-sufficiency, food production per capita, and power consumption of agricultural machinery per unit of cultivated land. Thus, these five indicators should be prioritized within the region to improve the synergistic safety of the WEF system.

WEF Stability
During the period 2010-2016, the stability of the WEF system in each prefecture-level city within the three northeastern provinces increased moving from south to north ( Figure  5). Over time, the stability of some prefecture-level cities in Jilin Province shifted to a safe state, with their energy self-sufficiency ratios ranging between 0.6 and 0. , specifically extensive energy consumption in the secondary industry in Jilin Province to revitalize the northeastern region, resulted in an increased rate of energy self-sufficiency in 2016. Most of the municipalities in Liaoning Province were in a state of insecurity relating to their stability, with their energy self-sufficiency rates ranging between 0.2 and 0.4, which were lower rates compared with those in the other two provinces. This result is linked to the amount of energy consumed in this large province, especially raw coal and crude oil, which account for 98% of the total energy consumption in the province.

WEF Stability
During the period 2010-2016, the stability of the WEF system in each prefecturelevel city within the three northeastern provinces increased moving from south to north ( Figure 5). Over time, the stability of some prefecture-level cities in Jilin Province shifted to a safe state, with their energy self-sufficiency ratios ranging between 0.6 and 0. with their energy self-sufficiency rates ranging between 0.2 and 0.4, which were lower rates compared with those in the other two provinces. This result is linked to the amount of energy consumed in this large province, especially raw coal and crude oil, which account for 98% of the total energy consumption in the province.

WEF Coordination
Contrasting with the spatial pattern of increasing stability of the WEF system moving from south to north across the three northeastern provinces (Figure 5), the coordination pattern was reversed, indicating an increasing trend moving from north to south ( Figure  6). Over time, there was an evident increase in the number of prefecture-level cities that demonstrated weak coordination in Heilongjiang and Jilin Provinces. The coordination levels of the WEF system in Baishan and Daxinganling indicated that they were extremely safe during the period 2010-2016, with the rate of agricultural water use in Baishan remaining consistently within a range of 0.06-0.1. In Daxinganling, the corresponding figures for agricultural water use were 0.76, 0.14, and 0.06 in 2010, 2013, and 2016, respectively, indicating continual improvements in the efficiency of water use in this area. In addition, national targets for the economic and social development of Heilongjiang Province in the 13th Five-Year Plan include an increase in the rate of water resources development and utilization up to 34% by 2020, and a resulting water supply capacity of 2 billion m 3 . These targets indicate that rational development and utilization of water resources should be accorded high priority.

WEF Coordination
Contrasting with the spatial pattern of increasing stability of the WEF system moving from south to north across the three northeastern provinces (Figure 5), the coordination pattern was reversed, indicating an increasing trend moving from north to south ( Figure 6). Over time, there was an evident increase in the number of prefecture-level cities that demonstrated weak coordination in Heilongjiang and Jilin Provinces. The coordination levels of the WEF system in Baishan and Daxinganling indicated that they were extremely safe during the period 2010-2016, with the rate of agricultural water use in Baishan remaining consistently within a range of 0.06-0.1. In Daxinganling, the corresponding figures for agricultural water use were 0.76, 0.14, and 0.06 in 2010, 2013, and 2016, respectively, indicating continual improvements in the efficiency of water use in this area. In addition, national targets for the economic and social development of Heilongjiang Province in the 13th Five-Year Plan include an increase in the rate of water resources development and utilization up to 34% by 2020, and a resulting water supply capacity of 2 billion m 3 . These targets indicate that rational development and utilization of water resources should be accorded high priority.
tively, indicating continual improvements in the efficiency of water use in this area. In addition, national targets for the economic and social development of Heilongjiang Province in the 13th Five-Year Plan include an increase in the rate of water resources development and utilization up to 34% by 2020, and a resulting water supply capacity of 2 billion m 3 . These targets indicate that rational development and utilization of water resources should be accorded high priority. Figure 6. Coordination of the WEF system in the three northeastern provinces (a-Daxinganling; b-Baishan). Figure 6. Coordination of the WEF system in the three northeastern provinces (a-Daxinganling; b-Baishan).

WEF Sustainability
The overall status of the sustainability of the WEF system in the three northeastern provinces was unsafe (Figure 7). Whereas the sustainability status of the WEF system in Tieling, Fuxin, Chaoyang, Huludao, and Dandong in Liaoning Province has remained very insecure over a period of several years, that of the city of Daqing in Heilongjiang Province has always been safer. The overall average GDP per capita of the three northeastern provinces was RMB 41,324, while the GDP per capita of Daqing remained within a range of RMB 81,325-130,707 from 2010 to 2016. The results of our comparative analysis revealed that a higher GDP per capita corresponded to a higher sustainability level of the WEF system and vice versa.

WEF Sustainability
The overall status of the sustainability of the WEF system in the three northeastern provinces was unsafe (Figure 7). Whereas the sustainability status of the WEF system in Tieling, Fuxin, Chaoyang, Huludao, and Dandong in Liaoning Province has remained very insecure over a period of several years, that of the city of Daqing in Heilongjiang Province has always been safer. The overall average GDP per capita of the three northeastern provinces was RMB 41,324, while the GDP per capita of Daqing remained within a range of RMB 81,325-130,707 from 2010 to 2016. The results of our comparative analysis revealed that a higher GDP per capita corresponded to a higher sustainability level of the WEF system and vice versa.

Levels of Synergistic Safety in WEFs
During the period 2010-2016, the synergistic safety index of the WEF system in most of the cities in the three northeastern provinces ranged between 0.4 and 0.6 ( Figure 8). This range of values indicated a state of boundary safety. The calculated weights of stability, coordination, and sustainability (Equation (5)) were 0.32, 0.36, and 0.32, respectively. These results suggest that the three WEF systems are mutually constraining and reinforcing, requiring equal attention. Furthermore, over time there has been a decline in the safety levels of these cities, especially in the northwestern part of the Songnen Plain. The security levels within the WEF system of cities within Liaoning Province indicated an un-

Levels of Synergistic Safety in WEFs
During the period 2010-2016, the synergistic safety index of the WEF system in most of the cities in the three northeastern provinces ranged between 0.4 and 0.6 ( Figure 8). This range of values indicated a state of boundary safety. The calculated weights of stability, coordination, and sustainability (Equation (5)) were 0.32, 0.36, and 0.32, respectively. These results suggest that the three WEF systems are mutually constraining and reinforcing, requiring equal attention. Furthermore, over time there has been a decline in the safety levels of these cities, especially in the northwestern part of the Songnen Plain. The security levels within the WEF system of cities within Liaoning Province indicated an unsafe state, with the security index values in nine cities remaining within a range of 0.2-0.4.

Analysis of the Degrees of Coupling and Coupling Coordination of the WEF System in the Three Northeastern Provinces
The spatial pattern of coupled coordination of the WEF system in the three northeastern provinces reflected an increasing trend moving from south to north ( Figure 9) and generally remained at a primary coordination. The coordination levels of prefecture-level cities in Heilongjiang and Jilin Provinces were in between the primary-intermediate levels in 2010 as well as 2013, and six prefecture-level cities in Heilongjiang were at the weakprimary level of coordination in 2016. In Liaoning Province, eight prefecture-level cities remained consistently at the weak-primary coordination level, and six prefecture-level cities were at the primary-intermediate coordination level from 2010 to 2016. The coupled coordination of the WEF system in the Tonghua shifted from barely primary to the primary-intermediate coordination level after 2013.

Analysis of the Degrees of Coupling and Coupling Coordination of the WEF System in the Three Northeastern Provinces
The spatial pattern of coupled coordination of the WEF system in the three northeastern provinces reflected an increasing trend moving from south to north ( Figure 9) and generally remained at a primary coordination. The coordination levels of prefecture-level cities in Heilongjiang and Jilin Provinces were in between the primary-intermediate levels in 2010 as well as 2013, and six prefecture-level cities in Heilongjiang were at the weakprimary level of coordination in 2016. In Liaoning Province, eight prefecture-level cities remained consistently at the weak-primary coordination level, and six prefecture-level cities were at the primary-intermediate coordination level from 2010 to 2016. The coupled coordination of the WEF system in the Tonghua shifted from barely primary to the primary-intermediate coordination level after 2013.
During the period 2010-2016, the coupling degree of the three northeastern provinces ranged between 0.8 and 1.0 (Figure 10), indicating a high level of coupling in which water, energy, and food strongly influenced each other. The coupling coordination degree fluctuated around a value of 0.6, which reflects a primary coordination level, indicating that the coupling coordination of the WEF system requires further improvement. In addition, we applied the GM (1, 1) model in conjunction with the coupling coordination (D) of the WEF system in the three northeastern provinces for the period 2010-2016 to predict D values for the period 2021-2026. The small probability error calculated using the posteriori difference test method was 0.86, and the variance ratio was 0.49, indicating that the GM (1, 1) model result was qualified. The development coefficient and the amount of the grey effect were 0.007 and 0.65 (Table 4). The D values ranged between 0.57 and 0.62 from 2017 to 2026, which reflected a weak coordination level, indicating that the coupling coordination of the WEF system requires further improvement (Table 5).
in 2010 as well as 2013, and six prefecture-level cities in Heilongjiang were at the weakprimary level of coordination in 2016. In Liaoning Province, eight prefecture-level cities remained consistently at the weak-primary coordination level, and six prefecture-level cities were at the primary-intermediate coordination level from 2010 to 2016. The coupled coordination of the WEF system in the Tonghua shifted from barely primary to the primary-intermediate coordination level after 2013. During the period 2010-2016, the coupling degree of the three northeastern provinces ranged between 0.8 and 1.0 (Figure 10), indicating a high level of coupling in which water, energy, and food strongly influenced each other. The coupling coordination degree fluctuated around a value of 0.6, which reflects a primary coordination level, indicating that the coupling coordination of the WEF system requires further improvement. In addition,  Table 4). The D values ranged between 0.57 and 0.62 from 2017 to 2026, which reflected a weak coordination level, indicating that the coupling coordination of the WEF system requires further improvement (Table 5).    In recent years, scholars have begun to explore the WEF system using different quantitative approaches. Niva et al. conducted spatial scenario analysis of China. The results manifested in the food sector, playing the leading role in the baseline water stress. The energy sector dominated the increases of the projected water stress index, and urbanization is projected to substantially affect the extent of water availability, especially in the eastern provinces [49]. Hua et al. calculated the water footprint of food and energy in 31 provinces in China to assess the consumption of water resources in food and energy production in different regions [62]. Their research mainly focused on macro-level evaluations and paid less attention to the external environment. In this study, we conducted a micro-level investigation focusing on the stability, coordination, and sustainability of 36 prefecture-level cities, and the results revealed a declining trend in the sustainability of the WEF system in the three northeastern provinces, with the level of synergistic safety fluctuating between 0.4 and 0.6, which is consistent with Li's research results [63]. Due to the different characteristics of the WEF symbiosis unit in different regions and the influence of the external environment on the WEF system being very different, the results are only appropriate for the three northeastern provinces, but the methods in this study can be applied to the global trends evaluation of China's economic sustainability. In addition, the Covid-19 pandemic may have some influence on the predictions for the next years in China. Among indicators, GDP per capita and energy consumption per 10,000 RMB of GDP may decrease. However, in this study we only considered internal flow of resources in 36 prefecture-level cities which were less affected by Covid-19 and did not involve import and export. Therefore, the trend of the results should show little deviation.

Suggestions on the Safety Development of the WEF System in the Three Northeast Provinces
High levels of water stress affect food security and induce increasing competition and potential conflict among sectors [14]. Globally, agriculture currently annually consumes 69% of annual water resources, and irrigated agriculture in particular is one of the main causes of water scarcity [64,65]. Therefore, improving the efficiency of water use in agriculture is an effective way of reducing water intensity [66]. For example, the cultivation of droughttolerant crops should be encouraged in Jilin and Heilongjiang Provinces. At the same time, agriculturalists should be made more aware of water consumption and conservation to enable them to optimize their water use [62], and the rational development and utilization of water resources should be prioritized in Heilongjiang Province. In addition, although the three northeastern provinces are richly endowed with energy resources, irrational exploitation of these resources and the adjustment of the industrial structure in recent years have led to mounting pressure on the region's energy supplies. The main factors influencing the safety of the WEF system in the northeastern provinces were agricultural water use efficiency, the energy consumption structure, and GDP per capita, with higher weight values indicating greater influence ( Figure 4). Consequently, the degree of dependence on external resources has gradually increased, causing a decline in the energy security evaluation index ( Figure 5). In particular, the energy consumption structure in Liaoning Province, which is inefficient and characterized by highly polluting energy consumption, needs to be transformed. Moreover, outdated production processes and methods and high energy-consuming equipment need to be eliminated and replaced. In addition, some regions should exploit their own resources and use them rationally, thereby increasing their GDP per capita.

Limitations
Comprehensive analysis of the safety development of the WEF system and the external social, economic, and natural factors may help to explain the relationship between WEF system and its external environment. However, there are still some limitations in this study. Firstly, when we collected data, in order to ensure the completeness of the data in 36 prefecture-level cities the type of data was limited, for example, in the food part we took into consideration only some food productions: cereals, potato, and soy, while we did not consider other factors (i.e., fruit, animal farming, and aquaculture as well as the sea fish industry). Secondly, the weight calculation method may also lead to the discrepancy compared with other research. Furthermore, due to the paucity of data for some prefecturelevel cities, the data of inputting the GM (1, 1) model were limited in this study, which also requires further refinement. Thirdly, the cross-departmental, multi-caliber, multi-scale nature of the WEF system and the ambiguity of system boundaries increase the difficulty of data acquisition and integration. Due to the limitations of data and quantitative methods, only the influence of some factors in the economy, society, environment, and land were considered in the evaluation process.
In future work, the "ecosystems," "biodiversity," "climate," and other relevant terms should be included in the database search. The comprehensive influencing factors and the dynamic feedback control simulation of the WEF system are also worth exploring.

Conclusions
In this study, which was premised on the symbiosis theory, we obtained a stability index, a coordination index, and a sustainability index for the WEF system. Using the comprehensive evaluation method, we constructed a synergistic safety index to evaluate the synergistic safety and coupling coordination of the WEF system in 36 prefecture-level cities in China's three northeastern provinces.
In the WEF system of the three northeastern provinces, the stability increased moving from south to north: Heilongjiang province was the most stable, while Liaoning province was the most unstable. The coordination had an increasing trend moving from north to south: Liaoning province had the highest coordination, while Heilongjiang province had the lowest coordination. The overall status of the sustainability in the three northeastern provinces was unsafe. The weights of the stability, coordination, and sustainability indexes for the synergistic safety of the WEF system were 0.32, 0.36, and 0.32, respectively, indicating equal emphasis is required to improve the integrated security status of the WEF nexus. The coupling degree ranged between 0.8 and 1.0, indicating a high level of coupling. The coupling coordination degree fluctuated around a value of 0.6, which reflected a primary coordination level, but it will range between 0.57 and 0.60 during the period 2021-2026, dropping to a weak coordinated level.  Institutional Review Board Statement: "Not applicable" for studies not involving humans or animals.
Informed Consent Statement: "Not applicable" for studies not involving humans.
Data Availability Statement: Data is contained within the supplementary material. The data presented in this study are available in supplementary material.

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