The Coupling Coordination Degree and Its Driving Factors for Water–Energy–Food Resources in the Yellow River Irrigation Area of Shandong Province

Abstract

Water resources, energy, and food are essential for the development of society, and they are strongly interdependent. The coupling and coordination relationships of the water–energy–food (WEF) system are important for regional resource security and high-quality development. The Yellow River Irrigation Area in Shandong Province, China, is a grain production base and has a substantial impact on national food security. To examine the water, energy, and food subsystem dynamics in this area, an evaluation system for the WEF system was established. A comprehensive weighting method based on game theory was employed to determine index weights. TOPSIS was used to assess the development level of the WEF system. A coupling coordination degree model was used to analyze the evolution of the coupling coordination degree of the WEF system from 2000 to 2020, and a GWR model was constructed to explore the spatial heterogeneity of its driving factors. The findings indicated that the development level of the WEF system in the study area was moderate, with a gradual upward trend. The coupling coordination degree fluctuated between 0.62 and 0.739. The GWR model revealed that temperature had an overall negative effect on the coupling coordination degree, with the greatest impact on the central irrigation area; the slope and NDVI had a negative effect, with increasing intensity from the southwest to the northeast; and rainfall had an overall positive effect, with the greatest impact on the irrigation area near the estuary in the northeast. Overall, the building area ratio had a negative effect on the coupling coordination degree, with exceptions in some areas. These research outcomes provide theoretical support for sustainable agricultural development in the Yellow River irrigation areas of Shandong Province and methodological reference data for studying collaborative resource utilization in irrigation regions.

1. Introduction

Water, energy, and food are essential resources for human survival and social development, and their security is closely related to sustainable human development [1]. With the increase in the global population, intensification of climate change, transformation of the natural environment, and continuous socioeconomic development, the three resource sectors of water, energy, and food are facing increasing demands and challenges. The growing demand for water, energy, and food is a concern for both developed and developing countries [2]. In the coming decades, the global demand for water, energy, and food will further increase. It is expected that by 2030, the demand for water resources, food, and energy will increase by 30%, 50%, and 40%, respectively [3]. As the world’s largest developing country, China faces certain risks and challenges for the future in terms of water, energy, and food in the face of international turbulence and global resource shortages. China’s water resources are scarce and unevenly distributed in time and space. China’s per capita water resources are approximately a quarter of the world average [4]. According to statistics from the Ministry of Water Resources of China, the annual agricultural production water shortage is approximately 26 billion cubic meters on the basis of the current typical water demand. According to the World Energy Statistical Yearbook (2019), China’s energy consumption has grown at an average rate of 3.9% over the past decade, accounting for 24% of global energy consumption. As the world’s largest energy-consuming country, China will still need a large amount of energy in the future, and the energy gap will continue to increase. According to the “China Rural Development Report (2021)” published by the Institute of Rural Development at the Chinese Academy of Social Sciences, by 2025, there will be a shortage of approximately 150 million tons of grain production in China that will need to be addressed. Considering the above problems, the coupling and coordination relationships of the WEF system should be explored from the perspective of WEF correlation to ensure the sustainable development of the region.

The relationships among water, energy, and food are strong and highly complex, and the production, processing, and transportation of food cannot be separated from water resources and energy [5]. This is because the production, storage, and consumption of energy consume water resources [6]; the extraction and scheduling of water resources consume energy [7]; and some grains can also be used as biomass energy [8,9,10]. Since the “Water-Energy-Food Security Linkage Conference” held in Bonn, Germany, in 2011, the water–energy–food nexus has become a popular research topic for politicians and scholars, and scholars worldwide have obtained many related results. Regarding WEF system evaluation, Mannan et al. [11] used the life cycle assessment method to analyze the levels of coordination between water resources, energy, and food; Mohammadour et al. [12] constructed a security evaluation index for the WEF system on the basis of the RAND pardee method and quantitatively evaluated the security of the WEF system; Li et al. [13] analyzed the risk levels of WEF systems on the basis of the copula function; and Li et al. [14] conducted a graded evaluation of the stability, coordination, sustainability, and collaborative security of WEF systems via a variable fuzzy evaluation method. In terms of model construction, Gu et al. [15] evaluated the safety of WEF systems on the basis of the pressure state response (PSR) model. Peng Shaoming et al. [16] proposed a layout plan for integrated optimization of water resource allocation, energy development, and grain production by constructing a collaborative optimization model for the WEF system. Hu et al. [17] used a system dynamics model to study the relationships among water, energy, and food in terms of urban metabolism. Wang Huimin et al. [18] combined the PSR model and system dynamics model to construct a model to explore the green development policy of the WEF system. Deng Peng et al. [19] constructed a coupled coordination degree model to study the evolution characteristics of the coupled coordination degree and promote regional sustainable development. At present, most scholars have evaluated the coupling and coordination of WEF systems in the time dimension; the lack of a hierarchy in research dimensions has made it difficult to thoroughly elucidate the coupling and coordination laws of WEF systems. Moreover, most of the existing research results are limited to regional scales, such as the national level [20], provincial level [21], and urban level [22], and few scholars have used irrigation areas as regional scales for research.

The Yellow River Irrigation Area in Shandong Province holds a crucial position in ensuring grain supply security for Shandong Province and ranks among the regions with the largest amount of water resources and agricultural energy utilization in China. Research on this area has predominantly focused on single-factor studies of water resources, including evaluations of irrigation water-saving techniques [23], water allocation planning and utilization [24], assessments of agricultural water security [25], and analyses of groundwater quality [26]. However, comprehensive analyses of the WEF system in this region are currently lacking. In the cropping systems of the Yellow River Irrigation Area in Shandong, crop rotation of wheat and maize has been widely adopted because of its significantly high yield. Nevertheless, the rapid increase in agricultural crop yields in the irrigation area has been accompanied by declining groundwater levels and substantial energy consumption. The advancing processes of urbanization and industrialization have led to further water, energy, and food security issues in the Yellow River Irrigation Area of Shandong, worsening the mismatch between the coordinated development of WEF subsystems and socioeconomic development. To address this situation, it is crucial to study the coupling coordination development patterns of the WEF system in the Yellow River Irrigation Area of Shandong Province. In this study, the Yellow River Irrigation Area in Shandong Province was selected as the research subject. A coupling coordination degree model and a GWR model were constructed to assess the coordination of water, energy, and food resource development over a period of years and investigate the spatial heterogeneity of their driving factors; the findings can contribute to the protection and sustainable development of agricultural resources in the irrigation area and the formulation of relevant policies. This effort was intended to provide theoretical references and insights for the scientific management and sustainable development of resources in the Yellow River Irrigation Area and Shandong Province as a whole.

2. Research Area and Data Sources

2.1. Research Area

The study area includes 28 Yellow River Irrigation Areas in Shandong Province, located at 34°33′–38°13′ N and 114°49′–119°06′ E (Figure 1). The total area is 388,300 km², currently including Liaocheng, Dongying, Binzhou, Jinan, Dezhou, Zibo, Tai’an, Jining, and Heze, and covers a total of 54 counties (districts) in 9 cities (prefectures) across Shandong Province.

The primary irrigation water source in the study area is the Yellow River, supplemented by small amounts of surface water from tributaries of the Yellow River, Hai River, and Huai River that serve the Shandong Yellow River Irrigation Area. This region experiences a warm, temperate monsoon climate characterized by mild temperatures and concentrated rainfall, with 60% to 70% of the annual precipitation occurring from June to August [27]. In recent years, the average rainfall in the lower reaches of the Yellow River has decreased due to climate change and human activities, leading to a significant reduction in natural river flow. The overall riverbed of the lower Yellow River shows signs of incision, and under same flow conditions, the water level of the Yellow River is generally decreasing. Additionally, water inflow in the upper reaches of the Yellow River has been insufficient, while the water demand for crop growth in the irrigation area continues to increase. Insufficient natural flow capacity in the irrigation area has worsened the increasingly prominent mismatch between the supply of and demand for agricultural irrigation water resources [25].

2.2. Data Sources

The research data primarily originated from authoritative sources such as the “Statistical Yearbook of Shandong Province”, the “Water Resources Bulletin of Shandong Province”, and statistical data from municipalities (counties, districts) within each irrigation area. In cases where data for specific years were missing, they were supplemented via weighted averaging or interpolation methods to ensure the completeness and accuracy of the study.

3. Method

This study presents a systematic framework for effectively exploring the coupling and coordination relationships among water, energy and food systems; the details are illustrated in Figure 2.

3.1. Comprehensive Evaluation Method

3.1.1. Construction of the Index System

The WEF system is a complex and interconnected mega-system with extensive content. The construction of an evaluation index system for assessing coupling coordination in this system involves consideration of uncertainties and diversity. By drawing on previous research [28,29,30,31,32], while comprehensively considering the interaction mechanisms among the water, energy, and food subsystems and strictly adhering to principles such as scientific rigor, systematicity, comprehensiveness, and data availability, a comprehensive evaluation index system for the WEF system in the Yellow River Irrigation Area of Shandong Province was developed. This index system comprises 18 evaluation indices (

Table 1. Comprehensive evaluation index system for the WEF system.

Target Layer	Subsystem	Index Layer	Index Attribute	Number
Comprehensive evaluation index system for the WEF system	Water resource system	Water resource utilization rate	Negative	A1
		Groundwater supply ratio	Negative	A2
		Agricultural water use ratio	Negative	A3
		Effective irrigation area ratio	Positive	A4
		Blue water use per unit area	Negative	A5
		Effective rainfall	Positive	A6
	Energy system	Fertilizer use per unit area	Negative	B7
		Electricity use per unit area	Negative	B8
		Plastic film use per unit area	Negative	B9
		Agricultural diesel use per unit area	Negative	B10
		Pesticide use per unit area	Negative	B11
	Food system	Wheat yield per unit area	Positive	C12
		Corn yield per unit area	Positive	C13
		Rice yield per unit area	Positive	C14
		Bean yield per unit area	Positive	C15
		Potato yield per unit area	Positive	C16
		Cotton yield per unit area	Positive	C17
		Yield of other grains per unit area	Positive	C18

).

Given the significant variations in cultivated land area across different irrigation areas in Shandong Province, average values across units of grain sowing area are considered in constructing the evaluation index system. This approach aims to reduce redundancy among evaluation indices and increase the accuracy of assessments, thereby ensuring the scientific validity and effectiveness of the comprehensive evaluation of the WEF system.

For the water resource subsystem, six indices closely related to water resource inputs were selected: the water resource utilization rate, groundwater supply ratio, agricultural water use ratio, effective irrigation area ratio, blue water use per unit area, and effective rainfall (A₁–A₆). The water resource utilization rate reflects the extent and potential of water resource development and utilization in the irrigation area. The agricultural water use ratio represents the current level of agricultural water use and the pressure it places on local water resources. The groundwater supply ratio reflects the degree of groundwater development and utilization. The effective irrigation area ratio is a key index for measuring the level of water infrastructure and agricultural production stability in the irrigation area. The blue water use per unit area reflects the utilization level of broad agricultural blue water resources in the irrigation area. The effective rainfall indicates the degree to which rainfall is directly or indirectly utilized during the growing season of cereal crops in the irrigation area, which is crucial for crop production planning and the determination of cultivation methods.

For the energy subsystem, this study focused on both direct (electricity, diesel) and indirect (fertilizers, pesticides, plastic film) inputs related to agricultural energy use in the irrigation area, and five indices were selected: fertilizer use per unit area, electricity use per unit area, plastic film use per unit area, agricultural diesel use per unit area, and pesticide use per unit area (B₇–B₁₁). These indices collectively reflect the intensity of different forms of energy inputs in agricultural production within the irrigation area.

For the food subsystem, seven indices closely related to grain production were selected on the basis of the actual grain production conditions in the irrigation area. These were the yields of wheat, corn, rice, beans, potatoes, cotton, and other grains per unit area (C₁₂–C₁₈). These indices collectively reflect the productivity of different types of crops within the irrigation area.

3.1.2. Calculation of the Weights

Subjective Weighting of the G1 Method

The G1 method, also known as the ordinal relationship analysis method [33], is a new subjective weighting method based on an improvement of the analytic hierarchy process (AHP) [34]. The G1 method assigns weights to different indices on the basis of an order relationship established by decision-makers. This method addresses the shortcomings of the AHP and reduces the computational effort of determining the weights of each index. It does not require the construction of a judgment matrix or a consistency check, making it easy to use. The advantages of the G1 method are its clear logic and simple calculation; the resulting weights are also sufficiently effective. This method uses three steps for determining subjective weights:

(1) Determine the order of importance of the evaluation indices:

$$c_1\succ c_2\succ \dots \succ c_n$$

(1)

where $n$ is the number of evaluation indices. The index on the left side of the symbol indicates greater importance.

(2) Determine the importance ratios between adjacent attributes in the order relationship:

$${\displaystyle \frac{w_{i-1}^1}{w_i^1}}=r_i(i=2,3,\dots ,\mathrm{n}),$$

(2)

The larger the $r_i$ value is, the greater the importance of the index $w_{i-1}^1$ compared with the index $w_i^1$. The reference values for $r_i$ are shown in

Table 2. Assignment of the relative importance between two neighboring indices.

ri	Meaning
1.0	Indices ci−1 and ci are equally important
1.2	Index ci−1 is slightly more important than index ci
1.4	Index ci−1 is clearly more important than index ci
1.6	Index ci−1 is much more important than index ci
2.0	Index ci−1 is extremely important compared with index ci

.

(3) Calculate the weights of each index:

$$w_i^1={(1+\sum _{i=2}^n\prod _{k=i}^n{r}_k)}^{-1},$$

(3)

$$w_{i-1}^1=w_i^1\cdot r_i(i=n,n-1,\dots ,2),$$

(4)

where $w_n^1$, $w_{i-1}^1$ and $w_i^1$ represent the weights of indices $c_n$, $c_{i-1}$, and $c_i$, respectively, obtained via the G1 method.

Objective Weighting via the Entropy Weighting Method

The entropy weighting method is an objective weighting method based on the theory of information entropy, which means that the smaller the information entropy and the greater the degree of dispersion, the more information is carried, and the greater the influence on the comprehensive evaluation [35,36,37].

The specific steps of applying the entropy weight method are as follows:

Assume that there are $m$ evaluation years and $n$ evaluation indices. The original data evaluation index matrix can be expressed as:

$$X=\left\{x_{ij}\right\}m\times n(j=1,2,3,\dots ,m;i=1,2,3,\dots ,n),$$

(5)

(1) Owing to the different dimensions of indices, it is necessary to standardize the data of each index:

$$Y_{ij}={\displaystyle \frac{x_{ij}-\min (x_{ij})}{\max (x_{ij})-\min (x_{ij})}},$$

(6)

$$Y_{ij}={\displaystyle \frac{\max (x_{ij})-x_{ij}}{\max (x_{ij})-\min (x_{ij})}},$$

(7)

(2) The characteristic proportions of each index are determined:

$$P_{ij}={\displaystyle \frac{Y_{ij}}{{\displaystyle \sum _{j=1}^m}{Y}_{ij}}},$$

(8)

(3) The information entropy value of the $i$-th index is solved:

$$E_i=-\ln {\left(m\right)}^{-1}\sum _{j=1}^m{P}_{ij}\times \ln P_{ij},$$

(9)

(4) The weights of each index are calculated via the following formula:

$$w_i^2={\displaystyle \frac{1-E_i}{n-{\displaystyle \sum _{i=1}^n}{E}_i}},$$

(10)

where $w_i^2$ represents the weight of index $x_i$ calculated using the entropy weighting method.

Comprehensive Weighting Based on Game Theory

The comprehensive weighting method based on game theory takes the Nash equilibrium as the coordination objective and aims to minimize the sum of the deviations between the combination weights and the basic weights; this makes the evaluation results more authentic and accurate [38,39,40].

(1) Build a basic weight vector set ${\omega }_k$:

If $G$ weight calculation methods are combined to assign weights to $n$ evaluation indices in the evaluation index system based on game theory, the corresponding weight vectors can be obtained as ${\omega }_k=\{{\omega }_{k1},{\omega }_{k2},\dots \dots {\omega }_{kn}\}(k=1,2,\dots ,G)$, and, furthermore, any linear combination of $n$ weight vectors can be obtained:

$$W=\sum _{k=1}^G{\alpha }_k{\omega }_k^T({\alpha }_k>0,k=1,2,\dots ,G),$$

(11)

where ${\alpha }_k$ is a weight coefficient.

(2) Build an optimized linear combination:

By optimizing the weight coefficients ${\alpha }_k$, the deviation between $\omega$ and ${\omega }_k$ can be reduced:

$$\min \Vert \sum _{k=1}^G{\alpha }_k{\omega }_k^T-{\omega }_i^T{\Vert }_2,$$

(12)

On the basis of the differential properties of the matrices, the first-order derivative conditions of the above optimization equation can be obtained:

$$\sum _{k=1}^G{\alpha }_k\cdot {\omega }_i\cdot {\omega }_k^T={\omega }_i\cdot {\omega }_k^T,$$

(13)

The linear equation system corresponding to the above equation is as follows:

$$\left[\begin{array}{cccc}{\omega }_1\cdot {\omega }_1^T& {\omega }_1\cdot {\omega }_2^T& \dots & {\omega }_1\cdot {\omega }_G^T\\ {\omega }_2\cdot {\omega }_1^T& {\omega }_2\cdot {\omega }_2^T& \dots & {\omega }_2\cdot {\omega }_G^T\\ ⋮& ⋮& \ddots & ⋮\\ {\omega }_G\cdot {\omega }_1^T& {\omega }_G\cdot {\omega }_2^T& \dots & {\omega }_G\cdot {\omega }_G^T\end{array}\right]\left[\begin{array}{c}{\alpha }_1\\ {\alpha }_2\\ ⋮\\ {\alpha }_G\end{array}\right]=\left[\begin{array}{c}{\omega }_1\cdot {\omega }_1^T\\ {\omega }_2\cdot {\omega }_1^T\\ ⋮\\ {\omega }_G\cdot {\omega }_1^T\end{array}\right],$$

(14)

The optimal linear combination $({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_G)$ obtained from the above equation is normalized to yield:

$${\alpha }^{\ast }={\displaystyle \frac{{\alpha }_k}{{\displaystyle \sum _{k=1}^G}{\alpha }_k}},$$

(15)

(3) Obtain the final combined weight value:

$$\omega =\sum _{k=1}^G{\alpha }^{\ast }\cdot {\omega }_k^T,k=1,2,\dots ,G,$$

(16)

3.1.3. TOPSIS Evaluation Method

The TOPSIS method, which is based on the distance between superior and inferior solutions, is a common comprehensive evaluation method. The basic principle is to find the best and worst solutions among limited options and then calculate the distances between each evaluation object and the best and worst solutions separately to obtain the relative closeness between each evaluation object and the optimal solution, which serves as the basis for evaluating the quality of the indices. This study used the TOPSIS method to evaluate the comprehensive development level of the WEF system and its subsystems in the Shandong Yellow River Irrigation Area. On the basis of the weights determined by the comprehensive game-theoretic weighting method, the TOPSIS method was used to rank the evaluation objects in terms of their strengths and weaknesses. The main calculation steps of this ranking method are as follows:

Assume that there are $n$ evaluation objects and $m$ standardized matrices for the evaluation indices:

$$\mathit{Z}=\left[\begin{array}{cccc}{z}_{11}& z_{12}& \dots & z_{1m}\\ z_{21}& z_{22}& \dots & z_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ z_{n1}& z_{n2}& \dots & z_{nm}\end{array}\right],$$

(17)

(1) Define the best and worst solutions for each index:

$${\mathit{Z}}^{+}=(z_1^{+},z_2^{+},\dots ,z_m^{+}), z_j^{+}=\max (z_{ij}),$$

(18)

$${\mathit{Z}}^{-}=(z_1^{-},z_2^{-},\dots ,z_m^{-}), z_j^{-}=\min (z_{ij}),$$

(19)

(2) Calculate the Euclidean distance from each evaluation object to the best and worst solutions as follows:

$$D_i^{+}=\sqrt{{\displaystyle \sum _{i=1}^n{\left({\omega }_j{z}_j^{+}-{\omega }_j{z}_{ij}\right)}^2}},$$

(20)

$$D_i^{-}=\sqrt{{\displaystyle \sum _{j=1}^m{\left({\omega }_j{z}_j^{-}-{\omega }_j{z}_{ij}\right)}^2}},$$

(21)

(3) Calculate the closeness of each evaluation to the optimal solution:

$$I_i={\displaystyle \frac{D_i^{-}}{D_i^{+}+D_i^{-}}},$$

(22)

where $I_i$ refers to the comprehensive evaluation index, which has a value range of [0, 1]. If the value is closer to 1, the result of the calculation is better, which here indicates a higher comprehensive development level of the WEF system or a subsystem.

3.2. Coupling Coordination Degree Model

System coupling, which means that two or more systems interact and influence each other, is a crucial concept in understanding complex systems. In this context, the coupling coordination degree model is adopted because it can comprehensively and synthetically consider the interaction relationship among the three systems of water resources, energy, and food. Moreover, it can clearly reflect the synergy of various elements within the system at different development stages, thus providing a significant basis for regional sustainable development planning. To explore the coupling coordination relationship between WEF systems, this work refers to existing studies [41,42,43,44], which have made significant contributions in analyzing the coupling relationship of multiple systems. Then, a coupling coordination degree model is established as follows:

(1) Calculate the coupling degree. The coupling degree refers to the degree of interaction and influence among various systems, and the formula is as follows:

$$C={\displaystyle \frac{3\sqrt[3]{I(w)\cdot I(e)\cdot I(f)}}{I(w)+I(e)+I(f)}},$$

(23)

In the above equation, $I(w)$, $I(e)$, and $I(f)$ refer to the comprehensive evaluation indices of the water resource system, energy system, and food system, respectively. $C$ is the coupling degree, where $0\le C\le 1$. The larger the value is, the greater the degree of correlation and interaction between the subsystems.

(2) Calculate the coupling coordination degree. This value reflects the degree of coupling coordination among various subsystems, indicating the quality of coordination. The relevant formulas are as follows:

$$T=\alpha I(w)+\beta I(e)+\gamma I(f),$$

(24)

$$D=\sqrt{C\times T},$$

(25)

In the above equation, $D$ is the coupling coordination degree; T is the comprehensive evaluation index of the three subsystems; and $\alpha$, $\beta$ and $\gamma$ are the weights of the water resource subsystem, energy subsystem and food subsystem, respectively. This study assumed that the impacts of the three subsystems on the development of irrigation areas were equally important, so $\alpha =\beta =\gamma =1/3$ was set. Referring to existing research [45,46,47], the specific level classification criteria for coupling coordination degrees are shown in

Table 3. Criteria for classifying the degree of coupling coordination.

Development Stage	Coupling Coordination Degree	Grade Standards
Extreme disorder	[0, 0.1)	Extreme dysregulation and decline
	[0.1, 0.2)	Severe dysregulation and decline
Basic disorder	[0.2, 0.3)	Moderate dysregulation and decline
	[0.3, 0.4]	Mild dysregulation and decline
Transition coordination	[0.4, 0.5)	On the brink of dysregulation and decline
	[0.5, 0.6)	Barely coordinated development
Moderate coordination	[0.6, 0.7)	Primary coordinated development
	[0.7, 0.8)	Intermediate coordinated development
High coordination	[0.8, 0.9)	Well-coordinated development
	[0.9, 1.0]	Highly coordinated development

.

3.3. Geographically Weighted Regression Model

The geographically weighted regression (GWR) model, possessing an extraordinary ability in spatial heterogeneity analysis, can efficiently capture the spatial variation characteristics of the WEF system’s performance in the study area, which are triggered by differences in natural conditions, economic development levels, and social factors. Simultaneously, it can thoroughly analyze the spatial disparities in the impacts of different factors on the WEF system, thereby providing precise information for formulating region-specific development strategies. Moreover, the GWR model incorporates the spatial and geographical location information of the data into the regression model. In contrast to the ordinary least squares (OLS) model, the GWR model provides more accurate local analysis results for the regression coefficients of each spatial unit and is better at explaining the influence of local areas in practical problems [48]. Therefore, this study adopts the GWR model to analyze the spatial heterogeneity of the driving factors of the WEF system.

In this study, the GWR model can be expressed as:

$$y_i={\beta }_0(u_i,v_i)+{\displaystyle \sum _{k=1}^p{\beta }_k(u_i,v_i)x_{ki}}+{\epsilon }_i (i=1,2,\dots ,n),$$

(26)

where $(u_i,v_i)$ represents the coordinate of the $i$-th basic research unit; ${\beta }_k(u_i,v_i)$ is the partial regression coefficient of the $k$-th key driving factor in the $i$-th basic research unit and is also a function of the geographic location of the sample point; $y_i$ is the coupling coordination degree of the WEF system in the $i$-th subregion of the study area; $x_{ki}$ represents the $k$-th key driving factor of the coupling coordination degree of the WEF system in the $i$-th region; ${\epsilon }_i$ is a random error term; and $p$ represents the number of explanatory variables.

According to the First Law of Geography, the observed values are weighted. When the geographic location $i$ changes, the assignment of the observation points also changes. The formula is as follows:

$$\mathit{\beta }(u_i,v_i)={\left[\begin{array}{c}{\mathit{X}}^T\mathit{W}(u_i,v_i)\mathit{X}\end{array}\right]}^{-1}{\mathit{X}}^T\mathit{W}(u_i,v_i)\mathit{Y},$$

(27)

where

$$\mathit{X}=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1z}\\ x_{21}& x_{22}& \cdots & x_{2z}\\ ⋮& ⋮& & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nz}\end{array}\right], \mathit{Y}=\left[\begin{array}{cccc}{y}_1& & & \\ & y_2& & \\ & & y_3& \\ & & & y_n\end{array}\right],$$

(28)

$$\mathit{W}\left(u_i,v_i\right)=\mathit{W}\left(i\right)=\left[\begin{array}{cccc}{w}_{i1}& 0& \cdots & 0\\ 0& w_{i2}& \cdots & 0\\ ⋮& ⋮& & ⋮\\ 0& 0& \cdots & w_{in}\end{array}\right],$$

(29)

$$\mathit{\beta }=\left[\begin{array}{c}{\beta }_0(u_1,{\nu }_1){\beta }_1(u_1,{\nu }_1)\cdots {\beta }_z(u_1,{\nu }_1)\\ {\beta }_0(u_2,{\nu }_2){\beta }_1(u_2,{\nu }_2)\cdots {\beta }_z(u_2,{\nu }_2)\\ \cdots \cdots \\ {\beta }_0(u_n,{\nu }_n){\beta }_1(u_n,{\nu }_n)\cdots {\beta }_z(u_n,{\nu }_n)\end{array}\right],$$

(30)

where $\mathit{\beta }$ is the estimator; $n$ denotes the number of basic spatial units; and $w_{in}$ is the weight of the $n$-th basic unit in geographic location $i$.

In the GWR model, using the Gaussian function as the spatial weight function yields better performance than using other functions. The Gaussian function is expressed as:

$$w_{ij}=\exp (-0.5{(d_{ij}/b)}^2),$$

(31)

where $d_{ij}$ represents the geographic distance between the basic research units $i$ and $j$; and $b$ is the bandwidth, which is an important parameter of the spatial weight function. The larger the bandwidth is, the smaller the magnitude of the decrease in the weight as the spatial distance increases.

In this study, the goodness of fit of the GWR model is measured at different bandwidths via the Akaike information criterion (AIC). The fitting effect of a model with a smaller $AIC$ value is better. If the $AIC$ difference between two models is >3, it indicates that there is a significant difference between the two models, and the model with a smaller $AIC$ value has a better fit. Assuming that the random error terms of the regression model follow an independent normal distribution, the $AIC$ expression of the regression model is as follows:

$$AIC=n\ln (RSS)+2q,$$

(32)

where $q$ represents the number of unknown parameters and where $RSS$ represents the residual sum of squares.

4. Result Analysis

4.1. Comprehensive Evaluation Index Analysis

4.1.1. Index Weight Results

On the basis of the methodology described in Section 3.1.2 of this article, the G1 method (sequential relationship analysis method) and entropy weight method were used to subjectively and objectively assign weights to 18 evaluation indices of the water resource system, energy subsystem, and food subsystem. According to Equations (1)–(4) and Equations (5)–(10), the subjective weight and objective weight of the WEF system can be obtained, respectively. A comprehensive weighting method based on game theory was subsequently applied to combine the subjective and objective weights obtained above. This method balances the advantages and disadvantages of different weighting methods to ensure that the combined weights are obtained scientifically and rationally. According to Equations (11)–(16), the comprehensive weight of the WEF system can be obtained. The final obtained combination weights are shown in Figure 3.

4.1.2. Analysis of the Comprehensive Development Level

According to Equations (17)–(22), the comprehensive development level of the WEF system in the study area can be calculated. Specific results can be seen in Figure 4 and Figure 5. The comprehensive development levels of the water resource subsystem, energy subsystem, grain subsystem, and WEF system in the Yellow River irrigation area of Shandong Province tended to increase to varying degrees and generally developed in a positive direction. Over time, the comprehensive development indices of the water resource subsystem, energy subsystem, and WEF system all exhibited significant and rapidly increasing indices. However, overall, the development levels of the water resource and grain subsystems lagged slightly behind that of the energy subsystem, especially the comprehensive development index of the grain subsystems, which grew relatively slowly. This inhibited the overall development rate of the WEF system to some extent.

Note: On the vertical axis, the areas from the Wangzhuang Irrigation Area to the Penglou Irrigation Area are the left-bank irrigation areas from the lower reaches of the Yellow River to the upper reaches in sequence; the areas from the Shuanghe Irrigation Area to the Yantan Irrigation Area are the right-bank irrigation areas from the lower reaches of the Yellow River to the upper reaches.

4.2. Analysis of Changes in the Coupling Coordination Degree

In this study, the comprehensive evaluation indexes of the three subsystems in the Yellow River diversion irrigation areas of Shandong Province obtained in Section 4.1 were combined with the coupling coordination degree model. Using Equations (17)–(25), the coupling coordination degree of the WEF system in the Yellow River diversion irrigation areas of Shandong Province was calculated. The calculated coupling coordination degrees of each irrigation area from 2000 to 2020 were obtained, and according to the classification standards in Table 3, these degrees were sorted into groups.

According to Figure 6 and the corresponding calculation results, from 2000 to 2020, the coupling coordination degree of the WEF system in the Yellow River Irrigation Area of Shandong Province fluctuated between 0.62 and 0.738, with an average value of 0.678, indicating moderate coordinated development; there was an overall upward trend. Regarding the spatial distribution, the coupling coordination degree of the irrigation areas on both sides of the lower reaches of the Yellow River in Shandong Province was greater than that of the irrigation areas on both sides of the upper reaches, which was consistent with the spatial distribution of the comprehensive evaluation index of the WEF system of the Yellow River diversion irrigation areas in Shandong Province. During the research period, the number of irrigation areas with primary coordination development continued to decrease, whereas the number of irrigation areas with intermediate coordination development continued to increase. In 2000, the coupling coordination degree of all 28 irrigation areas was of the primary coordination development type. However, by 2003, the coordination types of the nine irrigation areas on both upstream sides (the Yantan, Xiezhai, Liuzhuang, Susizhuang, Suge, Yangji, Chengai, Guonali, and Penglou irrigation areas) had decreased to barely coordinated development, whereas the remaining irrigation areas retained primary coordinated development. In 2008, the Bailongwan irrigation district was upgraded to intermediate coordinated development, whereas the other irrigation districts retained primary coordinated development. By 2013, the number of intermediate coordinated development-type irrigation areas on both sides of the downstream area had increased to nine (the Bailongwan, Bojili, Xingjiadu, Lijiaan, Panzhuang, Wangzhuang, Shuanghe, Mazhazi, and Liuchunjia irrigation areas), whereas the remaining irrigation areas retained primary coordinated development. By 2018, except for the nine irrigation areas on both upstream sides (the Yantan, Xiezhai, Liuzhuang, Susizhuang, Suge, Yangji, Chengai, Guonali, and Penglou irrigation areas), which were of a primary coordination development type, all irrigation areas were upgraded to an intermediate coordination development type. In 2019, except for the four irrigation areas of Suge, Chengai, Guo Nali, and Penglou, all irrigation areas exhibited intermediate coordinated development. The coordination type of the WEF system in the irrigation area in 2020 was the same as that in 2018. In summary, the coupling coordination degree of the nine irrigation areas on both upstream sides (the Yantan, Xiezhai, Liuzhuang, Susizhuang, Suge, Yangji, Chengai, Guonali and Penglou irrigation areas) fluctuated relatively strongly during the research period, but overall, it maintained an increasing trend. The coupling coordination degree of the other irrigation areas fluctuated less, and the development trend significantly improved.

4.3. Analysis of the Driving Factors of the Coupling Coordination Degree

In this work, when choosing the driving factors, 11 factors that posed potential risks to the WEF system proposed by the FAO in 2014 were considered, including population growth, urbanization rate, dietary diversification, cultural and social beliefs, climate change, governance, sectoral decision-making, international trade, industrial development, agricultural upgrading and technological innovation. Accounting for the distinctive geographical location and natural conditions of the Yellow River Diversion Irrigation Area in Shandong Province, by conducting a comprehensive review of multiple relevant materials and ensuring the representativeness of the selected driving factors and data accessibility, this paper selected five natural and social factors, namely, temperature, rainfall, slope, the normalized difference vegetation index (NDVI), and the building area ratio; their impact on the coupling coordination of the WEF system in the Yellow River Diversion Irrigation Area in Shandong Province was then explored. The specific details are shown in

Table 4. Driving factors of coupling coordination degree of the WEF system.

Variable	Dimension	Driving Factor	Number
Independent variable	Natural factor	Temperature	X1
		Rainfall	X2
		Slope	X3
		NDVI	X4
	Social factor	Building area ratio	X5
Dependent variable	Coupling coordination degree		Y

.

Since the GWR model examines the spatial heterogeneity of the driving factors, it is necessary to explore the spatial relationship of the coupling coordination degree of the WEF system before applying the GWR model. The spatial autocorrelation test was conducted via GeoDa (1.22.0.4) software, and the coupling coordination degree of the WEF system in the Yellow River Diversion Irrigation Area of Shandong Province showed clear spatial agglomeration. Considering the deficiencies in the fitting results of the OLS model and the significant spatial autocorrelation relationship of the coupling coordination degree of the WEF system in the Yellow River Diversion Irrigation Area of Shandong Province, the GWR model was employed to explore the spatial heterogeneity of the driving factors.

The comprehensive comparison results in

Table 5. Comparison of model results.

Model	AICc	R2	Adjusted R2
OLS	−137.4727	0.8298	0.7912
GWR	−141.4521	0.9562	0.9462

indicate that the AICc value of the GWR model is smaller than that of the OLS model, and the difference is 3.9794. When the numerical gap between AICc values is greater than 3, the fitting result of the model with the smaller value is considered better. Moreover, a comparison of the R² and adjusted R² values of the two models reveals that the results of the GWR model are superior, which again proves that the GWR model should be used to analyze the driving factors of the coupling coordination degree of the WEF system in the Yellow River Diversion Irrigation Area of Shandong Province.

The parameters of the GWR model were determined according to Equations (26)–(32), and the spatial heterogeneity of the driving factors was visualized via ArcGIS (10.8) analysis software. The final analysis results are shown in Figure 7, which represent the spatial heterogeneity effects of temperature, rainfall, slope, NDVI, and building area ratio on the coupling coordination degree of the WEF system in the Yellow River Diversion Irrigation Area of Shandong Province.

Figure 7a shows that the coefficient of the influence factor of temperature ranges from −0.077424 to 0.027784. In the vast majority of the irrigation districts in the study area, it is negative, indicating a negative correlation with the coupling coordination degree of the WEF system as a whole. Spatially, this coefficient shows a gradually decreasing trend from the southwest and northeast areas to the middle. Only the temperature influence factor coefficients of the Wangzhuang and Shuanghe irrigation districts in the northeast have a positive effect on the coupling coordination degree of the WEF system. The absolute value of the temperature influence factor coefficient of the irrigation district in the middle of the study area is relatively high, suggesting that the temperature has a stronger influence on the coupling coordination degree of the WEF system; that is, for each unit increase in temperature, the decrease in the coupling coordination degree of the WEF system is greater than those in the southwest and northeast irrigation districts.

Figure 7b shows that the coefficient of the driving factor of rainfall ranges from −0.000048 to 0.000219. In the vast majority of the irrigation areas in the study area, this coefficient is positive and has a positive correlation with the coupling coordination degree of the WEF system. That is, when other factors remain unchanged, the greater the rainfall, the higher the coupling coordination degree of the WEF. Only the Wangzhuang and Shuanghe irrigation areas in the northeast show a negative effect on the degree of coupling coordination of the WEF system. The coefficient of the driving factor of rainfall gradually increases from the southwest to the northeast.

Figure 7c shows that the coefficient of the driving factor of the slope ranges from −0.082329 to −0.000174. In all irrigation areas in the study area, this coefficient is negative, indicating a negative correlation with the coupling coordination degree of the WEF system globally. Spatially, this coefficient shows a gradually decreasing trend from the southwest to the northeast. The absolute value of this coefficient in the irrigation area near the estuary of the study area is relatively high, indicating that the influence of slope on the coupling coordination degree of the WEF system is relatively strong. That is, for each additional unit of slope increase, the decrease in the coupling coordination degree of the WEF system will be greater than that in the irrigation areas in the southwest and middle regions.

Figure 7d shows that the coefficient of the driving factor of NDVI ranges from −0.366856 to −0.052009. In all irrigation areas in the study area, this coefficient is negative, indicating a negative correlation with the coupling coordination degree of the WEF system globally. Spatially, it presents a gradually decreasing trend from the southwest to the northeast. The absolute value in the northeastern irrigation area of the study area is relatively high, suggesting that the influence of the NDVI on the coupling coordination degree of the WEF system is relatively strong. That is, for each additional unit of increase in the NDVI, the decrease in the coupling coordination degree of the WEF system will be greater than that in the central and southwestern irrigation areas.

Figure 7e shows that the coefficient of the driving factor of the building area ratio ranges from −0.043305 to 0.030691. In the vast majority of the irrigation areas in the study area, this coefficient is negative and has a negative correlation with the coupling coordination degree of the WEF system. The absolute value in the central irrigation areas is relatively high, indicating that in most of the central irrigation areas, a change in the building area ratio will cause a greater decline in the coupling degree of the WEF system. Only in a few irrigation areas, such as Penglou, Taochengpu, Guonali, Wangzhuang, and Shuanghe in the northeast and southwest, does the building area ratio have a positive effect on the coupling coordination degree of the WEF system.

5. Discussion

On the basis of the research results outlined in Section 4, it can be concluded that from 2000 to 2020, the comprehensive development level and coupling coordination degree of the WEF system in the 28 irrigation areas showed a stable upward trend. It is worth noting that both showed significant characteristics in which the downstream areas were superior to the upstream areas in space. In terms of the correlation between this pattern and the economic development status of the urban areas where the irrigation areas are located, the research results of this study are highly consistent with existing research results. Specifically, the lower reaches of the Shandong section of the Yellow River generally presented a relatively high level of economic development. The overall average economic growth indices of Dongying and Jinan were high, ranking first and second, respectively; Zibo and Binzhou followed closely behind. The economic development level in the upstream region was relatively low, and the overall average economic growth indices of Liaocheng and Heze cities ranked lowest, eighth and ninth, respectively [49]. In general, this result profoundly reveals that as important resources in human society, water, energy, and food have an inseparable relationship between their coupled and coordinated development and high-quality economic development as well as sustainable social development [47,50]. Specifically, the rational allocation and efficient utilization of water resources, energy resources, and food resources can not only promote the high-quality development of social economy but can also provide a solid material foundation for the sustainable development of society [51]. Therefore, in the process of promoting China’s economic and social development in the future, more attention should be paid to the coupled and coordinated development of resources to achieve the harmonious coexistence of economy, society, and environment [52].

The overall correlation between the temperature and the degree of coupling coordination of the WEF system is negative. The research findings of this paper closely align with existing research outcomes [53]. Temperature has a negative influence on water resources and agricultural output in the Yellow River Basin. The absolute value of the central irrigation area is relatively high. This is because the cultivated land area of the central irrigation area itself is relatively large. An increase in temperature will increase evaporation in the irrigation area, resulting in an increase in irrigation water demand. This may make already scarce water resources even scarcer, affecting the normal irrigation of crops and thereby threatening food production. Moreover, to meet the extraction, transportation and distribution requirements for irrigation water in irrigation areas, more energy is needed to drive equipment such as water pumps. Therefore, the temperature has a stronger effect on the coupling coordination degree of the WEF system in the central irrigation areas. Overall, rainfall has a positive correlation with the coupling coordination degree of the WEF system. The research results of this paper are consistent with existing results [54]. As rainfall gradually increases from inland areas to coastal areas, appropriate rainfall can increase surface runoff and groundwater resources in irrigation areas, provide more water for irrigation, promote the growth of crops, and improve the yield and quality of grain. Therefore, the positive impact of rainfall changes on the coupling coordination degree of the WEF is relatively significant. The slope has a globally negative correlation with the coupling coordination degree of the WEF system. The influence of slope on the coupling coordination degree of the WEF system in the northeastern irrigation area is much greater than that in the central and southwestern irrigation areas. This might be because the irrigation areas near the estuary are ecologically sensitive and the ecosystem is relatively fragile [55]. Minor changes in slope may trigger a series of ecological problems, such as wetland degradation and biodiversity reduction, which will indirectly affect agricultural production and water resource utilization. In addition, the soil in the irrigation areas near this area is mostly newly formed sedimentary soil, with fine particles, a loose structure and weak erosion resistance. An increase in slope will accelerate the soil erosion rate, and fertility will rapidly decrease, strongly affecting the growth of crops and leading to a significant decrease in the coupling coordination degree of the WEF system. The NDVI has a globally negative correlation with the coupling coordination degree of the WEF system. The influence of the NDVI on the coupling coordination degree of the WEF system in the northeastern irrigation area is much greater than that in the central and southwestern irrigation areas. This is possibly because in the northeastern irrigation area, an increase in the NDVI may mean a reduction in grassland and cultivated land areas and an increase in shrub and forest areas, resulting in a high demand for irrigation water; local water resources cannot meet this demand, leading to insufficient irrigation, a decline in grain production, and a greater impact on the coupling coordination degree of the WEF system [56]. Overall, the building area ratio has a negative correlation with the coupling coordination degree of the WEF system. This implies that when other factors remain unchanged, the larger the building area ratio is, the lower the coupling coordination degree of the WEF system. The influence of slope on the coupling coordination degree of the WEF system in the central irrigation area is much greater than that in the northeastern and southwestern irrigation areas. This might be because the industrial structure in the middle of the irrigation area of the Yellow River in Shandong Province is relatively simple and more sensitive to changes in the building area ratio. When the building area ratio changes, it may directly affect the layout and development of the core industry. In this case, each increase in the building area ratio has a significant effect on the coupling coordination degree of the WEF system.

In recent years, although the coordinated development of the WEF system in the Shandong Yellow River Irrigation Area has been continuously optimized, the comprehensive development of the water resource system and the grain system still lags behind that of the energy subsystem. Moreover, the coordinated development of WEF systems in various irrigation areas has also shown an uneven trend. Given the complex and interdependent relationship within the WEF system, the development status of any subsystem, whether it is progress or lagging behind, will have a profound impact on the coupled and coordinated development of the entire WEF system [57]. Therefore, to achieve the efficient, coordinated, and sustainable development of the WEF system in the Yellow River diversion irrigation areas of Shandong Province, it is necessary to comprehensively consider the internal connections among the three subsystems [58]. In the future development planning of irrigation districts, in order to ensure the sustainable development of the WEF system in this area, it is advisable that the widespread application of agricultural water-saving technologies in irrigation districts be actively and continuously promoted. Efforts should be dedicated to continuously improving the utilization efficiency of irrigation water in irrigation districts and accelerating the construction process of water-saving irrigation districts. Moreover, the planting structure of grain crops in irrigation districts should be continuously adjusted and optimized. The planting area of high-water-consuming crops should be moderately reduced, while the planting scale of high-yield and low-water-consuming grain crops should be increased. In addition, when formulating relevant policies, emphasis should be placed on considering the stability and safety of the energy system. Energy-saving agricultural machinery and equipment should be actively promoted. At the same time, efforts should be made to reduce energy consumption and energy utilization efficiency in agricultural production processes. Additionally, the application intensity of nitrogen fertilizers and pesticides should be reduced to reduce agricultural non-point source pollution. Communication and cooperation among water resource, energy and food management departments should be strengthened, and a sound cross-departmental consultation mechanism should be established. Finally, through the comprehensive use of various management means and technical strategies, involved parties should strive to consolidate and promote the coordinated, efficient and sustainable development of the WEF system in the Yellow River diversion irrigation districts in Shandong.

6. Conclusions

The coupling coordination degree of the WEF system is mostly concentrated at regional scales such as provinces, river basins, and the whole country. As a semi-artificial and semi-natural agricultural production system, the interaction relationships among water, energy, and food in irrigation districts are more complex and closer. Based on the conceptual framework of the WEF system’s nexus relationship, this study constructs an evaluation index system suitable for evaluating the coupling coordination level of the water–energy–food system in agricultural production in irrigation districts. Taking the Yellow River diversion irrigation districts in Shandong Province as an example, this study analyzes the coupling coordination degree and driving factors of the water–energy–food (WEF) system of 28 Yellow River diversion irrigation districts in the study area, reveals the evolution characteristics of the coupling coordination degree of the WEF system in the study area from 2000 to 2020, and explores the driving mechanism of the driving factors of the coupling coordination degree of the WEF system. The main conclusions are as follows:

(1) The comprehensive evaluation index of the WEF system in the Yellow River Irrigation Area of Shandong Province revealed that this area was at a moderate level of development overall and had a slow and stable upward trend during the research period. Specifically, both the energy subsystem and the water resource subsystem have made significant progress in their development, but compared with the energy subsystem, the development levels of the water resource and food subsystems have slightly lagged behind. In particular, the growth trend for the food subsystem was volatile, and the magnitude of change was relatively small, which, to some extent, indicated restricted overall development of the WEF system. In terms of spatial distribution, the irrigation areas on both downstream sides were generally better than the irrigation areas on both upstream sides in terms of the comprehensive development level of the WEF system, and the left-bank irrigation area also presented advantages over the right-bank irrigation area.

(2) The coupling coordination degree of the WEF system in the Yellow River irrigation area of Shandong Province fluctuated between 0.62 and 0.739, with an average value of 0.678, indicating a moderate coordinated development stage. In the long term, the overall development trend was stable and improving, with most irrigation areas achieving intermediate coordinated development. Regarding the spatial distribution, the coupling coordination degree of the WEF system in the 28 Shandong Yellow River Diversion Irrigation Areas showed clear regional characteristics, specifically, that both downstream-bank irrigation areas were superior to both upstream-bank irrigation areas.

(3) The results of the GWR model indicated that the temperature was generally negatively correlated with the coupling coordination degree of the WEF system, and the intensity of this effect increased from the northeast and southwest to the middle. The slope and NDVI were negatively correlated with the coupling coordination degree of the WEF system, and the intensity of this effect increased from the southwest to the northeast. Rainfall was generally positively correlated with the coupling coordination degree of the WEF system but was negatively correlated in a few irrigation areas in the southwest. Overall, the building area ratio was negatively correlated with the coupling coordination degree of the WEF system but was positively correlated with a few irrigation areas in the northeast and southwest.

These findings reveal the complex evolution characteristics of the coupling coordination degree of the WEF system in the Yellow River diversion irrigation areas of Shandong Province and the spatial heterogeneity of its driving factors, providing a scientific basis for the further optimization of irrigation resource allocation and development planning.