Assessing the Interaction Impacts of Multi-Scenario Land Use and Landscape Pattern on Water Ecosystem Services in the Greater Bay Area by Multi-Model Coupling

Abstract

Water ecosystem services (WESs) are intrinsically associated with the livelihood of urban residents and are frequently disrupted by human activities. Land use and landscape patterns are key driving factors of alterations in WESs. However, existing research primarily quantifies single-factor influences and often overlooks the interactions between these factors. This study addresses this gap by employing a multi-model coupling approach, integrating the Patch-generating Land Use Simulation (PLUS), Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model, and Geographical Detector (GD) models alongside various indicators to analyse the evolution of land use, landscape patterns and WESs in the Greater Bay Area from 2000 to 2020, and to simulate spatio-temporal change patterns in different scenarios from 2030 to 2050. Additionally, this study examines the multi-factorial interactions between land use, landscape patterns, and WESs. The results indicate that (1) urbanisation steadily increased, leading to intensified landscape fragmentation, and water yield (WY) and total phosphorus (TP) consistently increased, while total nitrogen (TN) in water gradually decreased; (2) urban areas exerted the most significant impact on WY in the Greater Bay Area while Patch density (PD) had a stronger influence on WY, and Shannon’s diversity index (SHDI) had the most pronounced effect on TN and TP; (3) the interaction between any two land-use types or landscape indices exerted a greater impact on WESs compared with the impact of individual factors alone. The interaction between urban areas and cropland substantially influenced WY ($\overline{q}$ = 0.634) and most strongly affected TN and TP in water ($\overline{q}$ = 0.74 and 0.73, respectively). SHDI and PD had the most significant impact on WY in the economic development scenario ($\overline{q}$ = 0.19) and exhibited the greatest influence on the TN and TP levels in the ecological priority scenario ($\overline{q}$ = 0.12 and 0.15, respectively). Our findings can provide theoretical and technical support for the integrated scientific planning of regional water ecosystems and the development of comprehensive land use policies in the future.

1. Introduction

Water ecosystem services (WESs) play a vital role in maintaining ecological balance, promoting ecosystem health and stability and providing essential services such as water yield (WY) and water purification (WP), which are indispensable for the sustainable development of human society [1,2]. However, the on-going processes of urbanisation and industrialisation have considerably impaired these services, posing a serious threat to human well-being and public health [3,4]. Existing studies indicate that water ecosystems are already experiencing significant degradation [5]. This degradation has resulted in a diminished capacity for WY and WP, thereby jeopardising the safety of drinking water supplies [3]. Given these challenges, the urgency of research on WESs is evident.

Over the past decades, the outward expansion of urban areas has directly led to a substantial reduction in water areas and severe degradation of the water ecosystem [6,7,8]. The contraction of water areas has diminished the capacity of the water ecosystem to supply essential services, severely impacting socio-economic and ecological environments [9,10]. In China, the continued expansion of urban areas and the consequent reduction in natural vegetation cover have led to decreased soil permeability, accelerated rainwater runoff and increased flood flows [11,12]. Consequently, research on the impacts of land-use change on WESs in urban agglomerations is essential. Landscape patterns considerably influence WESs, primarily owing to their role in controlling the flow of organisms and materials between landscapes [13,14]. Landscape fragmentation and connectivity have been further exacerbated by land-use changes, and these changes profoundly impact the level of WESs. The landscape patterns are often quantified via the landscape index method, and previous research has predominantly examined the effects of single factors on WESs using models such as geographically weighted regression, linear regression, and spatial econometrics [15,16,17]. For instance, Zhang et al. [17] explored the driving mechanisms involved in WESs in the Jianghuai region using the MGWR model, finding that SHDI exhibited significant positive and negative correlations with WESs, with PD exhibiting the strongest negative impact on WY. The interactions between ecosystem drivers and their impacts are also being investigated. Chen et al. [18] analysed the interactive effects of various land-use type pairs on the remote sensing-based desertification index (RSDI) in northern China from 2000 to 2018. Wu et al. [19] examined the influence of land cover on habitat quality in the Greater Bay Area (GBA) from 2000 to 2020, with a focus on the interactions between contributing factors. However, most of these studies focus on the single-factor effect or the interaction between factors in other ecological indicators, overlooking the comprehensive impact of interactions among multiple factors on WESs in the future.

Sustainable development of water ecosystems will face considerable challenges owing to national and social policies, highlighting the urgent need for modelling future WES scenarios [20,21]. Typically, these scenario simulations integrate land-use prediction models with ecohydrological models. Liu et al. [22] introduced the future land-use simulation (FLUS) model, while Liang et al. [23] developed the patch-scale land-use simulation (PLUS) model. The PLUS model employed the Random Forest Classification (RFC) algorithm to examine the relationship between multiple land-use expansions and multiple drivers. Meanwhile, the PLUS model can simulate more effective and realistic land use as well as landscape patterns by establishing a CA model based on multiple types of random patch seeds (CARS) [23,24]. In ecological modelling, the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model, developed by Stanford University, The Nature Conservancy and WWF, has been widely utilised. For instance, Wang et al. [25] used the PLUS-InVEST model to simulate the Han River Basin in 2050, revealing continued urban expansion and forest decline under a natural development scenario. The geodetector (GD) model is widely applied to analyse the influence of non-linear interactions between various factors on a target variable. Compared with previous studies, this study creates the PLUS-InVEST-GD research framework. This innovative framework is specifically designed to quantitatively evaluate the complex interactions between land use and landscape patterns, and their subsequent impact on WESs. By adopting this integrated approach, the study aims to provide deeper insights and more reliable data to inform the sustainable development of the WESs.

As one of China’s fastest-growing economic regions and a demonstration site for high-quality development, the GBA boasts the largest land area and population among the world’s major bay areas. WESs are critical to the economic development of urban agglomerations and the well-being of residents, particularly in terms of regional WY and WP. However, the land use and landscape patterns of urban agglomerations vary considerably over time and space during different development periods and, currently, limited research has been conducted on how their interactions will impact WESs in future scenarios. This study proposes a multi-model coupling approach to quantitatively assess the effects of multi-driver interactions on WESs. The PLUS, InVEST, and GD models were integrated with various indicators to analyse and simulate the spatio-temporal changes in land use, landscape patterns, and WESs in the GBA in different scenarios from 2000 to 2050. Subsequently, the complex mechanisms underlying the interactions between land use, landscape patterns and WESs were explored. The main contributions of our research include (1) a multi-model coupled approach created to study the mechanisms driving WESs by land use and landscape patterns; and (2) an exploration of the interactive effects of land use and landscape patterns on WESs. These findings significantly address existing gaps in the literature regarding WESs.

2. Materials and Methods

2.1. Study Area

The GBA is located in the south of China (111°21′–115°25′ E, 21°34′–24°23′ N), which comprises the nine cities of Guangzhou, Shenzhen, Zhuhai, Foshan, Huizhou, Dongguan, Zhongshan, Jiangmen, and Zhaoqing in Guangdong Province as well as the two special administrative regions of Hong Kong and Macao, with a total land area of 56,000 km² (Figure 1). The region experiences a subtropical humid monsoon climate, with an average annual temperature of 21–22 °C and annual precipitation ranges from 1600 to 2000 mm. By 2023, the economic output of the GBA surpassed 14 trillion yuan and its permanent population reached 86 million, making it one of the most open and economically vibrant regions in China. However, rapid urbanisation, continued expansion of urban areas and reduction in ecological land, including water, have considerably impacted the water ecosystem. Therefore, studying the WESs in the GBA is necessary and urgent.

2.2. Data Source

As can be seen from

Table 1. Data source in the study.

Data Type	Data Name	Resolution	Data Source
natural data	land use	1000 m	Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences
	soil type	1000 m
	precipitation	1000 m	WorldClimv2.1
	temperature	1000 m
	evapotranspiration	1000 m	National Tibetan Plateau Scientific Data Center
	DEM	30 m	Geospatial data cloud
socio-economic data	population	1000 m	Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences
	GDP	1000 m
	distance to primary, secondary and tertiary roads	30 m	National Center for Basic Geographic Information
	distance to the highway
	distance to the building
	distance to the railroad
	distance to the river

, the land-use data for 2000, 2010, and 2020 were obtained from the China Multi-Period Land Use Remote-Sensing Monitoring Dataset provided by the Chinese Academy of Sciences. This dataset primarily utilises Landsat remote-sensing images from the United States as its main information sources. The national-scale, multi-period land-use/land-cover thematic database of China was constructed through manual visual interpretation. To meet the research requirements, 25 secondary land-use classifications were reclassified and consolidated into six primary categories: cropland, forest, grassland, water, urban, and unused land. Additionally, we analysed 15 key drivers that considerably influence land-use change and WESs, including 10 socio-economic factors and 5 natural factors. The socio-economic factors include population, GDP and eight infrastructural characteristics (highways, railways, primary roads, secondary roads, tertiary roads, high-speed rail, buildings, and rivers), while the natural factors comprise temperature, precipitation, soil type, elevation and evapotranspiration. These driver data were resampled to new data with a resolution of 30 m. The data preprocessing work of this study is based on ArcGIS 10.7, including resampling and data registration.

2.3. Simulation Analysis of Land Use and Landscape Pattern

2.3.1. PLUS Model

The PLUS model is a cellular automaton (CA) framework designed to simulate patch-scale land-use/land-cover changes. This model integrates rule mining based on land expansion analysis and a multi-type random seed mechanism, enabling the identification of driving factors behind land expansion and the prediction of patch-level land use–landscape evolution. The PLUS model comprises two primary modules, namely the land expansion analysis strategy (LEAS) and the CA model with multi-class random patch seeds (CARS) [23].

The LEAS module extracts areas of land-use expansion between two historical periods and samples the expanded regions. Utilising the random forest classification algorithm, it identifies and quantifies the driving factors of various land-use expansions to determine the development probabilities of different land-use types [22,23,26].

The CARS module employs multi-variate random seed generation and a threshold reduction mechanism to simulate land-use changes based on the rules derived from the LEAS module, constrained by neighbourhood weights and a transfer matrix. Neighbourhood weight quantifies the difficulty of land-use conversion and is calculated using the following:

$${\mathsf{\Omega }}_{i,k}^t={\displaystyle \frac{con(c_i^{t-1}=k)}{n\times n-1}}\times w_k$$

(1)

where $W_k$ is the domain weight parameter, $n\times n$ represents the cell size, and $con(c_i^{t-1}=k)$ denotes the total number of grid cells occupied by land-use type kkk at the end of the cell iteration. Additionally, the domain weight of the instantaneous block type k at spatial unit i is expressed as ${\mathsf{\Omega }}_{i,k}^t$. The domain weight ranges from 0 to 1, where a higher value indicates greater expansion potential.

To analyse the future development trends of the GBA, this study integrates the current state of the region, economic development, and ecological red lines to construct three scenarios: the natural development scenario (NDS), the economic development scenario (EDS), and the ecological priority scenario (EPS). These scenarios simulate land-use changes in the GBA from 2030 to 2050. The scenarios are defined [19,27] as follows:

(1) The NDS assumes that land use in the GBA will follow natural development trends, meaning that, in the absence of significant policy interventions or external shocks, land use will continue to evolve according to historical patterns. Future land use under the NDS is predicted using Markov chain analysis based on historical land-use data from 2000 to 2010.

(2) The EDS reflects the transformation of land use due to rapid economic and urban growth in the region. In this scenario, the conversion of urban areas to other land uses is prohibited, while the probability of cropland, forest, grassland, water, and unused land to urban areas is increased by 20%.

(3) The EPS prioritises ecological protection in line with the Spatial Plan of Guangdong Province (2020–2035) and the ‘dual carbon’ target set in 2020. In this scenario, the probability of converting urban areas to cropland, forest, grassland, and water is increased by 30%, while conversions from these land uses to urban areas are prohibited.

2.3.2. Landscape Indices

In this study, landscape indices were calculated using FRAGSTATS 4.2, a widely recognised tool for landscape index analysis [28,29,30]. We selected 13 landscape indices (

Table 2. Landscape indices.

Type	Landscape Indices
Area-edge	Edge density (ED)
	Largest patch index (LPI)
	Total edge (TE)
Shape	Perimeter-area fractal represents dimension (PAFRAC)
Aggregation	Aggregation index (AI)
	Contagion index (CONTAG)
	Landscape shape index (LSI)
Subdivision	Number of patches (NP)
	Patch density (PD)
	Splitting Index (SPLIT)
	Landscape division index (DIVISION)
Diversity	Shannon’s diversity index (SHDI)
	Shannon’s evenness index (SHEI)

) that encompass area, shape, aggregation, and dispersion; these indicators comprehensively characterise landscape heterogeneity, complexity, structure, and fragmentation [30]. The moving window method was employed to calculate the landscape indices, which captures detailed characteristics and spatial variations in landscape patterns and allows for the analysis of temporal changes in the landscape [31]. After extensive testing, a window size of 1000 m × 1000 m was used to calculate the landscape pattern indices for the GBA.

2.4. WESs Assessment

WESs are categorised into two primary components: WY and WP. WY involves the quantitative estimation of runoff, calculated on a grid unit basis, and considers factors such as spatial variations in soil permeability and the evapotranspiration rates of different land-use types. WP refers to the capacity of the ecosystem to capture non-point source pollutants and mitigate water pollution through the functions of vegetation and soil. Total nitrogen (TN) and total phosphorus (TP) were chosen as indicators to assess the effectiveness of water quality purification. To quantitatively evaluate changes in WESs, the InVEST model was employed to compute these services for the GBA. In this paper, normal yearly precipitation and evapotranspiration during the study period has been mentioned in the Section 2.2, and the biophysical parameters used in this study were referred to related papers [32,33].

2.4.1. WY Capacity Assessment

WY is calculated using the WY module of the InVEST model, which is based on the Budyko curve and the principle of water balance [34]. The annual WY is determined using the following:

$$W{Y}_{xj}=\left(1-{\displaystyle \frac{AE{T}_{xj}}{PR{E}_x}}\right)\times PR{E}_x$$

(2)

where $W{Y}_{xj}$ represents the $WY$ for land-use type $j$, $AE{T}_{xj}$ denotes the actual evapotranspiration for land-use type $j$, and $PR{E}_x$ indicates the annual rainfall for the $x$ year.

2.4.2. WP Capacity Assessment

TN and TP are calculated using the nitrogen delivery ratio module of the InVEST model. This module estimates the TN to rivers based on land use and associated nitrogen and phosphorus load rates. TP is used as an indicator of water quality. The output of nitrogen and phosphorus is negatively correlated with the level of WP. The primary formula for these calculations is as follows:

$$AL{V}_x=HS{S}_X·po{l}_x$$

(3)

where $AL{V}_x$ represents the TN and TP of the pixel, $HS{S}_X$ is the nutrient output coefficient for nitrogen and phosphorus at the pixel level, and $po{l}_x$ denotes the hydrologic sensitivity fraction of the pixel $x$.

2.5. Geographic Detector

2.5.1. Factor Detection

The GD model is a statistical method used to identify spatial heterogeneity and its determinants, elucidating the impact of an independent variable X on a dependent variable Y. For this study, we employed factor detectors to analyse the landscape indices of water ecological services. The formula is as follows:

$$q=1-{\displaystyle \frac{1}{N{\sigma }^2}}{\sum }_{h=1}^L{N}_h{\sigma }_h^2$$

(4)

where $q$ represents the explanatory power of the landscape index, $h$ = 1, 2, 3, …, and $L$ denotes the classification or stratification of the landscape index. $N$ and ${\sigma }^2$ are the total sample size and variance, respectively, while $N_h$ and ${\sigma }_h^2$ represent the sample size and variance for the $h$ layer, respectively.

2.5.2. Interaction Detection

GD can evaluate the interaction between two different influence factors $X1$ and $X2$, determining whether their interaction enhances or diminishes the explanatory power of the dependent variable Y. This evaluation is typically performed by calculating the $q$ values for each factor individually, denoted as $q\left(X1\right)$ and $q\left(X2\right)$, and for their interactive combination, denoted as $q(X1\cap X2)$. By comparing these $q$ values, five types of interactions can be identified: (a) When $q(X1\cap X2)$ = $q\left(X1\right)+q\left(X2\right)$, the interaction is considered independent. (b) When $q(X1\cap X2)$ < $Min \left[q\right(X1), q(X2\left)\right]$, the interaction results in non-linear weakening. (c) When $Min \left[q\right(X1), q(X2\left)\right]$ < $q(X1\cap X2)$ < $Max \left[q\right(X1), q(X2\left)\right]$, the interaction shows single-factor non-linear weakening. (d) When $q(X1\cap X2)$ > $q\left(X1\right)+q\left(X2\right)$, the interaction leads to non-linear enhancement. (e) When $q(X1\cap X2)$ > $Max \left[q\right(X1), q(X2\left)\right]$, the interaction results in two-factor enhancement.

GD is a suite of statistical methods that detect spatial dissimilarities and reveal the underlying driving forces. GD is adept at analysing numerical data and identifying the interactive effects of two factors on a dependent variable. It uniquely discerns whether these factors have multiplicative interactions, a capability unmatched by other spatial interaction tools. Land use, landscape patterns, and WESs are all spatially quantified, providing an ideal context for applying GD to explore their interrelationships. In this study, we employed the GD method, facilitated by the geodetector data package in R, which primarily utilises the factor_detector and interaction_detector functions. The independent variables consist of categorical data of land use and landscape pattern indices, while the dependent variables encompass WY, TN, and TP.

3. Results

3.1. Changes in Land Use and Landscape Pattern During 2000–2050

3.1.1. Historical Land Use Analysis

The GBA underwent significant land-use transformation between 2000 and 2020 (Figure 2, Figure 3 and Figure 4), marked by a substantial increase in urban areas and a pronounced decline in forest, cropland and water. The period can be divided into two distinct phases: 2000–2010 and 2010–2020. During the first phase, rapid urban expansion and significant water loss were predominant, whereas other land-use types experienced relatively minor reductions. As shown in

Table 3. Land-use change from 2000 to 2020 (km²).

	Year	2000	2010	2020
Land Use Types
Cropland		14,439.94	12,632.43	12,092.92
Forest		30,607.81	30,021.74	29,674.57
Grassland		1221.88	1096.62	1182.57
Water		4384.05	4057.09	4012.72
Urban		4456.31	7314.19	8308.04
Unused land		23.43	11.36	6.54

, regions such as Jiangmen and Zhongshan experienced substantial water loss, primarily converted into urban areas and cropland. In the second phase, the pace of land-use change slowed. Urban areas continued to expand, albeit at a reduced rate of 99.39 km²·a⁻¹. Meanwhile, water decreased at a rate of 4.44 km²·a⁻¹, while grasslands began to recover at a rate of 8.6 km²·a⁻¹. The central part of the GBA, being more economically developed than the other regions, experienced the most significant urban expansion. Throughout the period, the total area of artificial surfaces increased by 2916.56 km², underscoring the accelerated urbanisation in the GBA. This expansion led to the continuous encroachment of artificial surfaces into water and other ecological land types, resulting in a significant loss of the natural environment.

3.1.2. Multi-Scenario Land-Use Change Modelling

Utilising land-use data between 2000 and 2010, the land use in 2020 was simulated using the LEAS and CARS modules within the PLUS model. The kappa coefficient of the simulation results when compared with actual land use in 2020 was 0.815, indicating that the PLUS model effectively captures spatial landscape changes and meets the requirements for predicting land use in this study. Consequently, based on the 2020 land-use data, land use was simulated in three scenarios for 2030–2050. The evolution of land use in these scenarios reveals distinct characteristics, as each scenario prioritises different environmental control objectives. The simulation results are presented in Figure 2.

The findings (Figure 3 and Figure 4) indicate that urban land expands in the NDS and EDS, leading to a reduction in water and other land-use types. Additionally, urban expansion is slower in the NDS than in the EDS, where the rapid development and expansion of urban areas disrupt the balance of land types. Conversely, in the EPS, water area increases for the first time, particularly in the western and north-eastern parts of the GBA, suggesting that the eco-green development policy enables better retention of water area.

3.1.3. Changes in Landscape Pattern at the Class Level

Table 4. Rate of change (%) in the landscape index at the class level from 2000 to 2020.

Land Use Type	Year	NP	PD	LPI	TE	ED	LSI	PAFRAC	DIVISION	SPLIT	AI
Forest	2000–2010	15.94	15.55	−0.96	−1.03	−1.37	−0.23	−1.49	0.06	1.57	−0.01
	2010–2020	−19.64	−19.57	−5.21	−1.44	−1.36	−0.63	2.28	0.38	9.44	0.00
Cropland	2000–2010	22.69	22.34	−24.68	−5.01	−5.34	1.48	−0.73	0.06	87.96	−0.43
	2010–2020	−18.24	−18.18	−33.61	−5.03	−4.95	−2.95	2.86	0.03	57.96	0.04
Grassland	2000–2010	4.53	4.10	−0.81	−5.32	−5.65	−0.27	−0.74	0.00	27.81	−0.46
	2010–2020	−4.80	−4.72	0.41	4.92	5.01	1.24	0.63	0.00	−14.60	0.22
Urban	2000–2010	−8.54	−8.85	503.5	20.98	20.57	−6.04	0.29	−0.12	−94.23	1.63
	2010–2020	−17.82	−17.77	−12.85	−10.49	−10.41	−14.98	−2.38	0.02	18.66	0.87
Water	2000–2010	3.81	3.48	−12.76	−2.27	−2.60	0.80	0.23	0.04	31.17	−0.24
	2010–2020	−20.70	−20.67	10.12	0.69	0.78	1.79	2.49	−0.02	−17.41	−0.22
Other	2000–2010	−10.77	−16.67	−50.00	−36.24	−36.58	−4.57	−2.31	0.00	310.52	−2.43
	2010–2020	−24.14	−20.00	−81.82	−19.92	−19.50	0.71	6.33	0.00	516.25	−3.44

demonstrates the changes in landscape indices at the class level for the period 2000–2020 in the GBA. In the first decade, forest patches exhibited the largest increase and decrease in NP (15.94%) and PAFRAC (−1.49%), respectively. Cropland patches exhibited the largest increases in SPLIT (87.96%), NP (22.69%), and PD (22.34%), with the largest decrease in LPI (−24.68%). Urban patches exhibited the largest increase in LPI (503.5%) and the largest decrease in SPLIT (−94.23%). Water patches showed the largest increase in SPLIT (31.17%) and the largest decrease in LPI (−12.76%). In the latter decade, forest patches demonstrated the largest increase and decrease in NP (19.64%) and SPLIT (−9.44%), respectively. Cropland patches exhibited the largest increase in SPLIT (57.96%) and the largest decreases in LPI (−33.61%), NP (−18.24%), and PD (−18.18%). Urban patches showed the largest increase in SPLIT (18.66%) and the largest decreases in NP (−17.82%) and PD (−17.77%). Water patches demonstrated the largest increase in LPI (10.12%) and the largest decreases in NP (−20.7%) and PD (−20.67%). Overall, the 2000–2010 period showed an increase in the agglomeration and less fragmentation of urban patches and an increase in the dispersion and fragmentation of forest, grassland, water, and cropland patches. The 2010–2020 period exhibited a decrease in the complexity and an increase in the dispersion of urban patches, a decrease in the fragmentation and an increase in the complexity of forest, grassland, and water patches and a decrease in the fragmentation of cropland patches, which were more discrete in their distribution.

It is demonstrated the changes in landscape indices at the class level in the GBA over the 2030–2050 period (

Table A1. Rate of change (%) in the landscape index at the class level from 2030 to 2050.

Land Use Type	Scenarios	Year	NP	PD	LPI	TE	ED	LSI	PAFRAC	DIVISION	SPLIT	AI
Forest	EDS	2030–2040	66.32	66.31	−0.55	3.51	3.51	4.35	1.96	0.05	1.58	−0.19
		2040–2050	20.91	20.91	−0.83	16.34	16.34	16.94	2.50	0.15	3.84	−0.67
		2030–2050	101.09	101.09	−1.38	20.43	20.43	22.02	4.51	0.20	5.47	−0.85
	EPS	2030–2040	−7.24	−7.24	1.22	4.42	4.42	3.48	0.73	−0.24	−4.19	−0.12
		2040–2050	−9.37	−9.37	0.76	0.05	0.05	−0.52	0.10	−1.76	−23.40	0.05
		2030–2050	−15.94	−15.94	1.98	4.47	4.47	2.95	0.83	−1.99	−26.61	−0.07
	NDS	2030–2040	130.61	130.60	0.79	91.19	91.19	91.02	10.71	−0.07	−1.65	−3.97
		2040–2050	10.87	10.87	5.35	7.60	7.60	8.41	0.90	−0.37	−8.33	−0.80
		2030–2050	155.67	155.68	6.18	105.72	105.72	107.09	11.71	−0.45	−9.84	−4.73
Cropland	EDS	2030–2040	27.08	27.07	−15.64	8.89	8.89	12.86	1.39	0.00	40.33	−2.45
		2040–2050	16.19	16.19	−1.94	−0.54	−0.54	2.27	0.66	0.01	20.79	−0.90
		2030–2050	47.65	47.65	−17.28	8.30	8.30	15.42	2.06	0.01	69.50	−3.33
	EPS	2030–2040	78.96	78.96	−7.25	5.27	5.26	8.93	0.99	0.01	35.57	−1.81
		2040–2050	31.76	31.75	−7.35	−1.11	−1.11	1.55	0.03	0.00	19.97	−0.70
		2030–2050	135.79	135.79	−14.06	4.09	4.09	10.62	1.02	0.01	62.64	−2.49
	NDS	2030–2040	75.21	75.21	−7.15	43.60	43.60	48.84	5.78	0.00	16.44	−10.06
		2040–2050	5.78	5.78	−11.12	5.41	5.41	8.41	0.73	0.01	26.11	−3.66
		2030–2050	85.35	85.35	−17.48	51.36	51.36	61.35	6.54	0.01	46.84	−13.35
Grassland	EDS	2030–2040	16.61	16.62	−10.81	2.61	2.61	5.55	1.68	0.00	16.34	−1.37
		2040–2050	8.66	8.66	−1.89	1.85	1.85	4.25	1.17	0.00	7.60	−1.17
		2030–2050	26.71	26.72	−12.50	4.51	4.51	10.03	2.87	0.00	25.18	−2.52
	EPS	2030–2040	10.93	10.92	6.87	5.44	5.45	4.21	0.44	0.00	−21.80	−0.99
		2040–2050	8.11	8.11	29.52	5.82	5.82	3.76	0.41	0.00	−17.58	−0.60
		2030–2050	19.93	19.92	38.41	11.58	11.58	8.13	0.85	0.00	−35.54	−1.59
	NDS	2030–2040	76.39	76.39	19.13	46.34	46.34	55.58	8.67	0.00	17.80	−19.14
		2040–2050	2.88	2.88	−31.82	−1.04	−1.04	2.34	0.50	0.00	39.59	−3.42
		2030–2050	81.47	81.46	−18.77	44.81	44.81	59.22	9.21	0.00	64.43	−21.90
Urban	EDS	2030–2040	8.04	8.04	19.28	13.76	13.76	7.49	1.53	−0.22	−28.89	−0.25
		2040–2050	3.03	3.03	6.01	11.61	11.61	7.06	1.20	−0.13	−14.34	−0.43
		2030–2050	11.31	11.31	26.46	26.97	26.97	15.07	2.74	−0.35	−39.09	−0.67
	EPS	2030–2040	−0.57	−0.55	0.00	0.01	0.01	0.01	0.50	0.00	0.00	0.00
		2040–2050	−0.41	−0.39	0.00	0.02	0.02	0.02	0.56	0.00	0.00	0.00
		2030–2050	−0.97	−0.94	0.00	0.03	0.03	0.03	1.06	0.00	0.00	0.00
	NDS	2030–2040	27.65	27.65	17.52	22.45	22.45	15.21	2.72	−0.26	−39.21	−1.18
		2040–2050	15.19	15.18	10.12	16.72	16.72	11.67	1.64	−0.12	−15.98	−1.05
		2030–2050	47.04	47.04	29.41	42.93	42.93	28.65	4.41	−0.38	−48.92	−2.21
Water	EDS	2030–2040	19.64	19.64	−40.11	1.60	1.60	4.21	1.74	0.05	133.05	−0.81
		2040–2050	8.03	8.02	−5.94	1.14	1.14	3.27	1.08	0.00	11.77	−0.68
		2030–2050	29.24	29.24	−43.66	2.76	2.76	7.62	2.84	0.05	160.48	−1.48
	EPS	2030–2040	6.21	6.21	2.26	5.14	5.14	3.27	0.65	0.00	−4.46	−0.19
		2040–2050	−7.76	−7.76	−0.72	0.28	0.28	−0.83	−0.14	0.00	1.06	0.25
		2030–2050	−2.03	−2.03	1.52	5.43	5.43	2.42	0.51	0.00	−3.45	0.05
	NDS	2030–2040	44.37	44.37	1.38	17.34	17.34	19.02	4.02	−0.01	−2.27	−3.98
		2040–2050	7.60	7.60	29.33	1.19	1.19	2.99	0.66	−0.02	−32.97	−1.16
		2030–2050	−2.03	−2.03	1.52	5.43	5.43	2.42	0.51	0.00	−3.45	0.05
Other	EDS	2030–2040	55.06	51.72	0.00	−5.08	−4.97	8.45	3.52	0.00	64.67	−10.99
		2040–2050	25.71	27.27	0.00	−14.69	−14.71	−6.44	−1.18	0.00	15.98	−1.76
		2030–2050	94.94	93.10	0.00	−19.02	−18.94	1.47	2.30	0.00	90.99	−12.56
	EPS	2030–2040	46.02	43.75	0.00	−10.93	−11.05	2.30	1.28	0.00	91.40	−9.52
		2040–2050	12.06	13.04	0.00	−11.38	−11.44	−3.63	−0.67	0.00	21.82	−2.30
		2030–2050	63.64	62.50	0.00	−21.06	−21.22	−1.41	0.60	0.00	133.16	−11.60
	NDS	2030–2040	52.66	51.61	0.00	−6.82	−6.92	8.44	3.55	0.00	78.83	−12.30
		2040–2050	43.41	42.55	0.00	−7.66	−7.77	0.22	0.36	0.00	6.72	−5.00
		2030–2050	118.93	116.13	0.00	−13.96	−14.15	8.68	3.92	0.00	90.84	−16.69

). In the EDS, forest and grassland patches show the largest increase in PD (101.09% and 26.71%, respectively) and the largest decrease in LPI (−1.38% and −12.5%, respectively). Cropland and water patches exhibit the largest increase is in SPLIT (69.5% and 160.48%, respectively) and the largest decrease is in LPI (−17.28% and −43.66%, respectively). Urban patches show the largest increase in ED (26.97%) and the largest decrease in SPLIT (−39.09%). In the EPS, cropland patches demonstrate the largest increase in PD (135.79%) and the largest decrease in LPI (−14.06%). Grassland patches show the largest increase in LPI (38.41%) and the largest decrease in SPLIT (−35.54%). Water patches exhibit the largest increases in ED (5.43%) and the largest decreases in SPLIT (3.45%). Forest patches exhibit the largest decreases in SPLIT (−26.61%) and PD (−15.94%). In the NDS, forest, cropland, grassland, water, and urban patches show the largest increases in NP and PD and a continuous decrease in SPLIT and AI. Comparing the three scenarios, the EDS and NDS indicate an increase in patch fragmentation, a decrease in agglomeration and an increase in complexity for all land-use types. Conversely, the EPS shows a decrease in fragmentation and an increase in agglomeration for forest and urban patches and an initial increase in fragmentation followed by an increase in agglomeration for water patches.

3.1.4. Changes in Landscape Pattern at the Landscape Level

The most significant increases at the landscape level in the GBA from 2000 to 2010 were observed in NP (5.65%) and PD (5.29%), while the most significant decrease was observed in CONTAG (by 1.88%) (

Table A2. Rate of change (%) in the landscape index at the landscape level from 2000 to 2020.

Landscape Indices	2000–2010	2010–2020	2000–2020
NP	5.65	−17.60	−12.95
PD	5.29	−17.53	−13.17
LPI	−0.96	−5.21	−6.12
TE	0.01	−3.55	−3.55
ED	−0.34	−3.47	−3.80
LSI	−0.20	−3.28	−3.47
PAFRAC	−0.67	1.84	1.17
CONTAG	−1.88	−0.41	−2.28
DIVISION	0.05	0.40	0.45
SPLIT	1.10	9.19	10.39
SHDI	3.51	1.15	4.70
SHEI	3.51	1.15	4.70
AI	0.01	0.12	0.13

). During the 2010–2020 period, the most substantial increase was recorded in SPLIT at 9.19%, while the most significant decreases were observed in NP and PD, declining by 17.6% and 17.53%, respectively. Other indices exhibited reductions in the following order: LPI (−5.21%), TE (−3.55%), ED (−3.47%), and LSI (−3.28%). Landscape diversity and fragmentation increased, while connectivity and complexity declined from 2000 to 2010. Conversely, landscape diversity and dispersion increased, accompanied by reductions in landscape fragmentation from 2010 to 2020.

As shown in

Table 5. Rate of change (%) in the landscape index at the landscape level from 2030 to 2050.

Landscape Indices	EDS			EPS			NDS
	2030–2040	2040–2050	2030–2050	2030–2040	2040–2050	2030–2050	2030–2040	2040–2050	2030–2050
NP	17.32	8.24	26.99	13.52	5.03	19.23	61.43	9.78	77.22
PD	17.32	8.24	26.99	13.52	5.03	19.23	61.43	9.78	77.22
LPI	−0.55	−0.83	−1.38	1.22	0.76	1.98	0.79	5.35	6.18
TE	8.70	7.92	17.31	4.56	0.08	4.64	45.94	8.04	57.68
ED	8.70	7.92	17.31	4.56	0.08	4.64	45.94	8.04	57.68
LSI	8.53	7.78	16.98	4.45	0.07	4.53	45.17	7.95	56.71
PAFRAC	1.63	1.38	3.03	1.03	0.21	1.25	5.77	0.94	6.76
CONTAG	−0.75	−0.90	−1.64	−0.03	0.25	0.22	−6.26	−1.37	−7.55
DIVISION	−0.10	0.01	−0.09	−0.24	−1.75	−1.99	−0.35	−0.54	−0.88
SPLIT	−2.37	0.27	−2.10	−3.99	−22.33	−25.43	−6.84	−9.74	−15.91
SHDI	−0.25	−0.33	−0.58	−0.63	−0.42	−1.04	−0.24	−0.25	−0.49
SHEI	−0.26	−0.32	−0.58	−0.63	−0.42	−1.05	−0.25	−0.25	−0.50
AI	−0.73	−0.73	−1.46	−0.32	−0.01	−0.32	−4.52	−1.21	−5.67

, the most significant increase at the landscape level in the GBA in the EDS from 2030 to 2050 occurred in PD, with a rise of 26.99%, and a slight decrease in SPLIT by 2.1%. In the EPS, the largest increase at the landscape level was observed in NP and PD, with an increase of 19.23%, accompanied by a significant decrease in SPLIT of 25.43%. In the NDS, NP, and PD showed the most substantial increase (by 77.22%), while SPLIT demonstrated the largest decrease (by 15.91%). Overall, all three scenarios demonstrate an increase in patch fragmentation and complexity, an increase in aggregation and a decrease in diversity at the landscape level.

3.1.5. Analysis of Spatial and Temporal Variations in Landscape Indices Under a Moving Window

Our testing revealed that a window size of 1000 m yields a more stable landscape index and accurately reflects the spatial information patterns. Five landscape indices were selected to analyse the temporal and spatial changes in the landscape pattern of the GBA (Figure A1 and Figure A2). CONTAG, SHDI, and PD exhibited an increase followed by a decrease from 2000 to 2020, with high values primarily concentrated in the central and southern regions. Conversely, LSI, AI, and ED showed modest increases, with high values predominantly observed in Zhaoqing in the north-west, as well as in southern Guangzhou and northern Huizhou. The central core areas of Guangzhou, Foshan and Dongguan within Shenzhen generally showed decreases in AI and LSI and increases in CONTAG, ED and PD from 2030 to 2050. In 2050, the EPS exhibited higher values of AI and SHDI and lower values of PD and ED compared to the other two scenarios. By contrast, the NDS and EDS showed higher values of CONTAG and PD. Overall, the EPS demonstrated increased aggregation and decreased fragmentation in the north-west and eastern regions, while the EDS and NDS exhibited increased fragmentation, higher complexity and reduced diversity in the central and southern regions.

3.2. Changes in WESs from 2000 to 2050

Table 6. Statistical table of changes in WESs from 2000 to 2050.

WESs	2000	2010	2020	2050 (NDS)	2050 (EDS)	2050 (EPS)
WY(×108 m3)	529.4	606.9	530.5	579.8	584.5	541.9
WN(t)	27,939.05	27,011.02	26,942.56	24,250.72	25,563.06	21,541.79
WP(t)	2200.42	2280.44	2350.38	2390.39	2537.63	1920.83

and Figure 5 illustrate the temporal and spatial variations in WY and WP in the GBA from 2000 to 2050. WY increased from 529.4 × 10⁸ to 606.9 × 10⁸ m³ between 2000 and 2010, marked by an expansion of high-value areas in the central region, a reduction in high-value areas in the south-eastern region. TN decreased from 27,939 to 27,011 t, with a reduction in TN in high-value areas in Guangzhou, Dongguan and Shenzhen. TP increased from 2200.42 to 2280.44 t, with low-value areas expanding in Jiangmen and Zhuhai and high-value areas growing in Foshan and Zhongshan. WY decreased from 606.9 × 10⁸ to 530.5 × 10⁸ m³ between 2010 and 2020, with a notable increase in low-value areas, particularly in the western region of the GBA. TN decreased from 27,939 to 26,942.56 t, with a decline in TN in high-value areas across the GBA, while TP increased from 2280.44 to 2350.38 t, with an increase in TP in high-value areas, particularly in the eastern region. In the 2050 NDS, WY increases from 530.5 × 10⁸ to 579.8 × 10⁸ m³, with the increase in WY in high-value areas primarily concentrated in the central region, while TP increases from 2350.38 to 2390.39 t, with high-value areas exhibiting TP increase in Shenzhen and the south-western region. In the 2050 EDS, WY increases from 530.5 × 10⁸ to 584.5 × 10⁸ m³, with high-value areas exhibiting WY increase primarily in Jiangmen and the central part of the GBA. TN decreases from 26,942.56 to 25,563.01 t, with a reduction in TN in high-value areas across the GBA, particularly in Guangzhou. TP increases from 2350.38 to 2537.63 t, with a decrease in TP in low-value areas, particularly in Huizhou. In the 2050 EPS, WY increases from 530.5 × 10⁸ to 541.9 × 10⁸ m³, with a decrease in WY in low-value areas, particularly in Jiangmen. TN decreases from 26,942.56 to 21,541.79 t, with high-value areas exhibiting diminishing TN in Guangzhou and Dongguan, while low-value areas exhibiting increased TN in Zhaoqing. TP decreases from 2350.38 to 1920.83 t, with an overall reduction in high-value areas across the GBA.

3.3. Analysis of Single Factor Detection Results

3.3.1. Single-Factor Detection Results of Land Use Impact on WESs

A t-test with a significance level of 0.05 was conducted for WY and WP (Figure 6). The results indicate that urban, forest, water, and cropland patches had a significant effect on the spatial distribution of WY from 2000 to 2050 (p < 0.05). The effect of urban patches on WY was significantly higher than those of the other land-use types, with a multi-year average q-value of 0.48. The q-value first increased and then decreased from 2000 to 2020, reaching approximately 0.5 or higher in 2050. Forest patches also considerably affected the spatial distribution of WY, particularly in the NDS and EDS, where the q-mean value exceeded 0.3. By contrast, the q-mean value for forest patches was 0.2 in the EPS. The influence of water patches on WY ranged from 0.1 to 0.2, peaking at 0.25 in 2020. The q-values for cropland and other land types were all less than 0.1, suggesting that water patches were depleted in these areas due to evaporation or other processes. The effect of grassland patches on WY was minimal in the EPS and insignificant in the EDS and NDS. At the same time, the results demonstrate that all land-use types significantly affected the spatial distribution of TN and TP (p < 0.05). Cropland (with a multi-year average q-value of 0.63) and forest patches (with a q-value of 0.4) were the land-use types with the greatest effect on TN. Forest patches had the highest influence on TP, with a multi-year average q-value of 0.5.

3.3.2. Single-Factor Detection Results of Landscape Pattern Impact on WESs

Figure 7 illustrates the impact of different landscape pattern indices on WY and WP from 2000 to 2050. First, a t-test at a significance level of 0.05 was conducted for WY. The results reveal that the five landscape indices significantly affected the spatial distribution of WY in 2000, 2010, and 2050 (p < 0.05) but not in 2020 (p > 0.05). Generally, the effects of landscape indices on WY were modest, with the exception of the 2050 scenarios (EDS and EPS), where the q-values for PD, LSI, and ED were notably high, approaching 0.04 to 0.1. Second, a t-test at a significance level of 0.05 was applied to WP. The results indicate that all landscape indices significantly affected the spatial distributions of TN and TP (p < 0.05). The influence of each landscape index on WP was generally below 0.1. SHDI and CONTAG had the most significant impact on TN, with multi-year average p values of 0.05 and 0.04, respectively. SHDI also had the greatest impact on TP, with a multi-year average p value of 0.04.

3.4. Analysis of Interaction Detection Results

3.4.1. Interaction Detection Results of Land Use Impact on WESs

The results of the interaction test of GD between land-use types and WY for the period 2000–2050 (Figure 8) indicate that the interaction effects of each pair of land-use types on WY improved relative to the single land-use type and exhibited a somewhat regular pattern over time. The q-values for interactions between land-use types from 2000 to 2050 were all below 0.5. Among these interactions, eight groups demonstrated two-way enhancement, while seven groups exhibited non-linear enhancement. Between 2000 and 2010, the interaction between urban patches and cropland showed the highest q-values, reaching 0.33 and 0.43, respectively. In 2020, the interaction between urban and water patches exhibited the highest q-value, reaching 0.34. In 2050, the q-values for interactions among the three scenarios showed significant enhancement. In the EDS, the q-values for interactions between urban patches and the other five land-use types exceeded 0.6, with seven groups showing non-linear enhancement. Notably, the interactions between urban patches and cropland showed the highest q-value (0.68). In the EPS, the q-values for interactions between urban patches and the other five land-use types were all above 0.5. Nine groups demonstrated dual-factor enhancement, and six groups showed non-linear enhancement. Particularly, the interactions between urban, cropland and water patches exhibited higher q-values. In the NDS, the q-values for interactions between urban patches and the other five land-use types were above 0.6. Eight groups showed dual-factor enhancement, and seven groups exhibited non-linear enhancement. Furthermore, the interactions between urban patches and cropland showed the highest q-value (0.69). These results indicate that interactions between urban patches and other land-use types significantly increase the effect on WY.

The results of the interaction test of GD between land-use types and WP for the period 2000–2050 (Figure 9 and Figure 10) indicate that the interaction effects of each pair of land-use types on WP improved compared to single land-use types and exhibited a somewhat regular pattern over time. For TN, the interaction q-values between cropland and the other five land-use types exceeded 0.5, with the highest values observed for cropland interacting with urban areas, reaching 0.77, 0.76, 0.66, 0.76, 0.76, and 0.73. The interaction q-values for forest and water exceeded 0.5, with values of 0.62, 0.58, 0.50, 0.52, 0.58, 0.52, and 0.58 for forest and 0.52, 0.53, 0.59, and 0.55 for water. These results indicate that the effect on TN is significantly enhanced when cropland interacts with other land-use types. For TP, the dual-factor interaction q-values of forest and cropland with other land-use types generally decreased over time. However, three pairs of interactions consistently showed high q-values: cropland and urban areas, forest and water and forest and cropland. The q-values for cropland and urban interactions and for forest and water interactions were all above 0.6, while those for forest and cropland interactions were above 0.5. These results indicate that the effect on TP is significantly greater when cropland interacts with urban and when forest interacts with other land-use types.

3.4.2. Interaction Detection Results of Landscape Pattern Impact on WESs

Figure 11 illustrates the results of the interaction test between the five aforementioned landscape pattern indices and the GD of WY from 2000 to 2050. The results demonstrate that the interaction effect between pairs of landscape indices on WY was enhanced compared to single landscape indices, exhibiting a somewhat consistent interaction pattern across different periods. The interaction q-values of the landscape indices exhibited only minor enhancements, all remaining below 0.024, with 10 groups showing a twofold enhancement from 2000 to 2020. Between 2000 and 2010, the largest interaction q-values were observed for SHDI and ED, reaching 0.022. In 2050, the interaction q-values were significantly enhanced across the three scenarios. In the EDS, interactions between SHDI and ED and between SHDI and LSI exhibited markedly high q-values of 0.213 and 0.212, respectively. Three groups showed dual-factor enhancement, seven groups exhibited non-linear enhancement, and seven interactions had q-values of ≥0.15. In the EPS, the q-values for interactions between SHDI and PD, ED and LSI were 0.129, 0.123 and 0.123, respectively. These interactions included one group with dual-factor enhancement, nine groups with non-linear enhancement, and six interactions with q-values exceeding 0.1. In the NDS, interaction results were less enhanced than the previous scenarios. One group showed dual-factor enhancement, and nine groups exhibited non-linear enhancement. The highest q-value was 0.021 for the interaction between SHDI and ED, while the q-values for interactions between SHDI and LSI and between PD and ED reached 0.019. The results indicate that the impact of SHDI on WY was significantly enhanced when interacting with PD, ED, and LSI, with the effect being most pronounced in the EDS.

Figure 12 and Figure 13 present the results of the interaction tests between the five landscape pattern indices and WP from 2000 to 2050. The results indicate that the interaction effects of each pair of landscape pattern indices on WP were enhanced compared to individual indices and demonstrated a consistent pattern in different periods. For TN, SHDI exhibited significant interaction effects with the other four landscape indices, particularly with CONTAG, PD, and ED, showing the highest q-values of 0.044, 0.052, 0.045, 0.087, 0.118, and 0.066. Notably, in the EPS, the interaction q-values of SHDI with the other four landscape indices all exceeded 0.1, indicating a substantial improvement. These results demonstrate that the interaction of SHDI with other landscape indices significantly enhances its effect on TN and WP, with the most pronounced increase observed in the EPS. SHDI and CONTAG consistently showed higher interaction q-values with the other indices for TP from 2000 to 2050. In the 2050 EDS and EPS, five interactions had q-values exceeding 0.1. These results indicate that interactions involving SHDI and CONTAG significantly enhance the effect on WP, particularly in the EPS and EDS.

4. Discussion

4.1. Land Use and Landscape Interactions on WES

Although previous studies have explored the impacts of changes in land use and landscape patterns on WESs, most studies have been unifactorial, lacking the detection of interactions between factors, particularly in future scenarios. Given that changes in WESs are influenced by multiple factors [35,36], this study applied GD to identify and analyse these interactions. The findings reveal that the effects of interactions between different land-use types on WESs in the GBA were generally more significant than the effects of individual indicators alone. Our results align more closely with those of Liang et al. [33], suggesting a similar trend in the data that interactions between urban and other land-use types had a substantial impact on changes in WY. The most notable effects were observed when urban areas interacted with cropland and water, resulting in improved WY owing to the conversion of significant areas of cropland and water into urban areas. Furthermore, interactions between cropland, forest, and other land-use types significantly influenced the spatial distribution of TN and TP. Interactions with urban areas, in particular, had the most pronounced effect, increasing TN and TP, which may contribute to the deterioration of water quality [37]. Additionally, the interaction of the landscape pattern index, CONTAG, with PD, and ED had the most significant impact on WY, highlighting that landscape connectivity and the degree of patch fragmentation and complexity are critical for regional WY [38]. The interaction between the SHDI and CONTAG had the most significant effect on WP. The diversity and connectivity of the landscape promoted the export flow of N and P, influencing the region’s WP [39]. In contrast, our findings are inconsistent with the results of Zhang et al. [40], which may be attributed to different research scales and variable definitions, that in urban areas with rapid urbanisation, the impact of urbanisation indicators on WESs was not always significant while, in rural areas, urbanisation indicators had a significantly negative impact. These differences highlight the complexities within the field and suggest the need for further investigation to reconcile these divergent outcomes.

4.2. Effect of Land Use and Landscape on WESs in Multiple Scenarios

Our findings align with several previous studies that report significant effects of land use and landscape pattern changes on WY [41,42,43] and WP [44,45]. Urban expansion in the GBA was found to significantly influence changes in WY, particularly in the future NDS and EDS. In these scenarios, large areas of cropland and forest were converted to urban areas, resulting in a dramatic increase in hardened surfaces, which had the most pronounced effect on WY. This increase in WY is primarily due to the hardened surfaces preventing water infiltration, leading to more direct surface runoff [42]. By contrast, cropland, forest, and grassland are restored in the EPS and the expansion of urban areas slows, resulting in a lower overall impact on WY. The restoration of forest and grassland increases soil moisture retention, leading to more significant changes in WY [46]. Cropland and forest had the most substantial effects on TN and TP, particularly in the EPS and NDS. An increase in cropland exhibited a stronger influence on water quality, aligning with studies that suggest that increases in cropland and urban areas lead to decreased WP [47]. Moreover, changes in landscape patterns were found to significantly impact WESs in different periods, although such changes did not fundamentally alter dominant relationships [17]. The most significant effects of PD, LSI, and ED on WY were observed in 2050. PD reflects landscape fragmentation and heterogeneity, which play a critical role in WY [14,38]. This indicates that the complexity of patch shapes also strongly influences regional WY. SHDI had the most substantial impact on TN and TP, and it reflects landscape diversity and heterogeneity, with patches characterised by higher SHDI composition reducing the impact of pollutants on water quality [39]. This suggests that landscape pattern diversity plays a crucial role in the degradation of TN and TP. CONTAG significantly affects WY, TN, and TP, as it responds to patch connectivity. Enhanced connectivity promotes material transfer and ecological processes, thereby increasing WY and TP [48].

4.3. Policies and Implications

The development policy in the EPS fosters favourable conditions for the restoration of water. However, as urbanisation progresses, the trend of urban expansion remains irreversible, posing a significant threat to the sustainable development of the region [24,49,50]. Currently, the GBA faces challenges related to high rates of water resource development and utilisation, regional water scarcity and widespread water quality issues. Based on the findings of this study, the following policy recommendations are proposed:

(1) The GBA should strictly adhere to the three control lines—ecological protection red lines, basic cropland, and urban development boundaries—in the national land-use spatial planning. This will protect water and forest from degradation, strictly control the scale of urban expansion and reserve space for future urban development, thereby mitigating the impact of urban expansion on WESs [51,52]. In rapidly urbanising areas such as Guangzhou, Dongguan, and Foshan, WESs are particularly vulnerable, with WY and WP capacity at risk of deterioration.

(2) The blue-green spatial pattern [53] should be optimised through the development of green ecological infrastructure [54] within urban agglomerations. Efforts should be made to prevent landscape fragmentation and reduce shape complexity while maintaining the diversity and heterogeneity of landscapes. These measures will enhance WY, improve WP capacity and strengthen regional WESs.

4.4. Limitations

Furthermore, this study acknowledges certain limitations inherent in the modelling approaches employed. For instance, the PLUS model’s assumptions regarding the probabilities of future land use transitions and the GD model’s determination of the sample size both exert a degree of influence on the experimental outcomes. These considerations necessitate further contemplation in subsequent research endeavours to bolster the validity and reliability of the findings.

5. Conclusions

As a crucial component of the socio-economic-natural composite ecosystem [55,56,57], water ecosystems play a vital role in identifying the factors affecting their service functions and their interactions [47,48]. However, the interaction mechanism of land use and landscape patterns on WESs remains unclear. In this study, a novel multi-model coupling method was proposed, integrating the PLUS, InVEST, and GD models with various indicators, to quantitatively assess the effects of historical and future multi-scenario land use and landscape pattern interactions on WESs in the GBA from 2000 to 2050. The findings are summarised as follows: (1) The GBA has been undergoing a rapid transformation of land-use type, such as water and cropland, into urban areas, accompanied by increasing patch fragmentation and decreasing diversity. WY showed a trend of initial increase, while TN exhibited a steady annual decline and TP increased annually. (2) Urban areas have the greatest impact on WY, while cropland and forest have the most significant influence on TN and TP. Moreover, interactions between different land-use types exceed the effects of individual factors. Furthermore, the interaction between urban areas and water has a pronounced effect on WY, whereas the interaction between cropland and forest is most significant for TN and TP. (3) PD, ED, and LSI exerted the most considerable influence on WY, while SHDI and CONTAG most significantly affected WP. Notably, the interaction between SHDI and PD had the greatest effect on WY, while the interaction between SHDI and CONTAG had the most significant impact on WP. Overall, through the study of the interactions between land use, landscape patterns, and WESs in different scenarios, we evaluated the hydrological and ecological effects of landscape patterns and land use. This study elucidates the underlying response mechanisms of WESs to shifts in landscape patterns and land use, offering a more comprehensive perspective for enhancing future WES levels. It is anticipated to significantly influence future land-use planning and landscape pattern development in the GBA, serving as a valuable scientific reference.