Research on the Change in Prediction of Water Production in Urban Agglomerations on the Northern Slopes of the Tianshan Mountains Based on the InVEST–PLUS Model

Abstract

Assessing how land use change will affect water production ecosystem services is essential to developing sound water resource management and ecosystem conservation. The results of a coordination analysis of land-use intensity and water yield based on future land-use simulation projections are useful for future land-use planning. To effectively assess water production rates, the PLUS and InVEST models were used to dynamically assess the changes in water production occurring in the urban agglomeration on the northern slopes of the Tianshan Mountains from 2000 to 2030 under different scenarios of land-use change. The results show that the water-production rates in the study area from 2000 to 2020 were 517.26 × 10⁶ m³, 582.28 × 10⁶ m³, and 456 × 10⁶ m³, showing an increasing and then decreasing trend, with the water production function decreasing from the foothills of the Tianshan Mountains to the north and south, with values of 509.10 × 10⁶ m³³, 510.90 × 10⁶ m³, and 502.28 × 10⁶ m³ being presented for the three scenarios in 2030. The rapid economic development scenario presents the lowest water yield values and the ecological conservation development scenario presents the highest water yield values. Changes in water production rates are closely related to changes in land use, which can be verified further by an analysis of the coordination between land- use intensity and water production. For this study area, the ecological conservation development scenario may be more in line with the future urban development pattern, and the results of the present study can provide some scientific references for land-use planning.

1. Introduction

Water is the basis for human survival and social development, and the relationship between people and water has an important impact on the coordination of human–earth systems [1]. With the increase in water demand and water shortages occurring due to the deterioration of the water environment, there is a threat to the water-holding capacity of regional ecosystems [2], which severely hinders food production and energy development [3]. In order to achieve sustainable development goals and maintain harmony between human beings and land systems, the trade-off between urban development and water ecosystem security is crucial [4,5]. As an important way in which ecosystems participate in hydrological regulation processes, water yield refers to the process and capacity of ecosystems to store and hold water in a given time and space [6], the drainage that occurs during the non-rainy season, and the recharge of the ecological base flow of a watershed to ensure water for human life and social development [7]. Quantitative assessments and visualizations of water yield are therefore increasingly observed as means of elucidating human–water interactions [8], which can help us to achieve the effective management and conservation of water resources, with important implications for sustainable development and human well-being at a regional scale [9].

Water production reflects the renewable capacity of regional water resources and the carrying capacity of the water environment. The contradiction between water supply and demand poses an increasing challenge to human production and life [10]. Regional water-retention capacity is related to vegetation, species survival, and various aspects of social life, and it is necessary to explore the delineation of water-saving ecological space under the guidance of water resources management in order to alleviate the conflicts between water use, water environments, water ecology, and resource development [11,12]. In recent years, many scholars have conducted a large number of studies at different levels regarding the relationship and mechanisms between human activities and water resource changes [13,14]. Water resource changes at the basin scale are significantly influenced by climate change and human activities [15]. At the urban scale, the impact of human activity on the water-producing functions of the ecosystem is more evident, and its effects are mainly reflected in the different land-use scenarios that change the formation of the watershed. Land-cover types change according to people’s needs, and these alter the physical and chemical properties of the substrate, ultimately affecting hydrological components, such as surface runoff, soil-water content, and evapotranspiration, and threatening the water security of the area [16,17]. Models are commonly used to assess water production, and the SWAT, PLOAD, HSPF, and InVEST models are commonly used in research practice [18]. Taking into account the geographical characteristics of the study area, the InVEST model is more applicable to the water production simulation of this study than former three models [19]. In addition, the model takes into account natural conditions, human activities [20], land-use [21], and water use trade-offs [22], making it more comprehensive. The InVEST model has been used to assess water production in geographically diverse areas, such as Hawaii, Indonesia, and China [23,24]. In previous studies, it has been observed that changes in water yield values by site type over recent years in the study area are generally consistent with the trends in area change [25], suggesting that climate and land-use changes are the main reasons for the decline in water yield services [26]. This also suggests that spatial variations in land-use driven by human activities are an important factor influencing the spatial and temporal variations in water yield. In practice, parameter settings need to be corrected in order to integrate them into the decision-making process. An important parameter that must be entered into the InVEST model when estimating water yield is the Z-parameter. Z is an empirical constant that reflects local precipitation patterns and hydrogeological characteristics, with typical values ranging from 1 to 30. The calibration of the Z-parameter can also be performed by comparing model data with observed data. Its primary role in the model is to influence estimates of water yield based on local precipitation patterns and hydrogeological characteristics. However, it is determined by comparing observed streamflow data, which can lead to a high degree of subjectivity and is unlikely to reflect the surface characteristics of the catchment [27]. Water production is an important ecosystem service that provides the basis for human survival and development, and water production distributions must be calculated and mapped with scientific accuracy. We have done extensive experiments in setting Z-parameters. Comparisons were made with government water resource bulletins to ensure the accuracy of the parameters. Predicting changes in the ecosystem of water production provides new ideas for studying future land-use change and the evolution of water production. Land-use change provides a clearer understanding of the ecological changes and their characteristics resulting from the interaction between human activities and the natural environment. The impact of land-use change on ecosystem services can be better reflected using land-use scenario simulations. In 2021, a raster-based Patch Generation Land-use Simulation (PLUS) model was developed to reveal the potential drivers and their different contributions to change. Rules for land use change were identified to improve the model’s ability to simulate realistic landscape patterns. The model has become one of the common support methods used by policymakers and managers. Scholars have used the model to conduct research [28]. These studies provide a valuable reference for regional water resources management and ecological planning. We found that, in addition to considering the impact of land type conversion on changes in water yield, considering also the harmonious relationship between land-use intensity and water yield functions in arid and semi-arid regions would facilitate a more systematic and comprehensive assessment of the impact of land-use change on regional water yield. In addition, the scientific validity of the model output depends on the accuracy of the input parameters. This study proposes a feasible and replicable framework for land use decision-making in ecologically fragile areas based on land use modeling techniques and a coordination degree model. It not only contributes to an in-depth understanding of the future urban growth patterns of the Tianshan North Slope urban agglomeration but also helps to inform different land-use stakeholders. The study will reveal the underlying patterns of water production evolution in the Tianshan north slope urban agglomeration, driven by the “Belt and Road” construction, to provide scientific reference for the sustainable development of the Tianshan north slope urban agglomeration and for the environmental protection of its ecosystem.

The overall objective of this study is to project future land-use changes under three different scenarios in 2030 and to analyze the impact of land-use on the water-producing functions of urban agglomerations to achieve a balance between economic growth and ecological conservation in ecologically fragile areas. The specific objectives of the study are (1) to project future land-use and water production changes in 2030 under three scenarios and assess the mechanisms of their interaction by coupling the InVEST and PLUS models, and (2) to determine the combined impacts and relative importance of land-use changes on changes in ecosystem services.

2. Overview of the Study Area

The urban agglomeration on the northern slopes of the Tianshan Mountains is located in the southern area of the Junggar Basin, between 42°78′ and 45°59′ N, 84°33′ and 90°32′ E, south of the Tianshan Mountains and north of the Gurbantunggut Desert (Figure 1; the mapping space referenced is WGS_1984_UTM_Zone_45N). The main natural vegetation consists of grassland and shrubs, and the woodland is dominated by evergreen coniferous forests in the southern mountains. The main soil types are brown calcareous soils, wind-sanded soils, salt soils, irrigated silts, and gray desert soils, and glacial and permafrost soils are widely distributed in the middle- and high-mountain areas. Under the combined influence of human activities and the natural climate, the salt desertification of the oasis margins in the region is significant, glacier retreat is severe, and the fragile ecological environment in the region can constrain sustainable socio-economic development.

3. Data Sources and Research Methodology

3.1. Data Sources

The data covered in this study include land-use data from 2000 to 2020 for the formation of Tianshan North Slope. The drivers of land-use change include natural conditions, accessibility, and socio-economic factors, of which elevation, slope and slope direction, average annual temperature, and rainfall data were selected for the natural conditions. Accessibility factors include distance vector data from rivers, railways, highways, national roads, provincial roads, county roads, etc., obtained by calculating Euclidean distances. Socio-economic factors were selected as the population density and nighttime light index, which were extracted using NPP/VIIRS nighttime light data (from the NOAA website). NDBSI data for the urban agglomeration on the northern slopes of the Tianshan Mountains were calculated based on MODIS Aqua surface reflectance data using MYD09A1, and NDVI data were calculated using the Google Remote Sensing cloud computing platform, calling the annual average data calculated from the 2020 continuity index MYD13A1V6 a 500-m resolution product. The above data were transformed into 1 km × 1 km raster data after a series of data pre-processing work such as projection transformation, cropping, and resampling in the ArcGIS analysis software.

3.2. Research Methodology

The framework of this study is presented in Figure 2. Firstly, historical LUCC data were used to predict future land use scenarios through the PLUS model, followed by the estimation of water yield in the study area through the InVEST model, and finally the analysis of the factors influencing the change in water yield values.

3.2.1. Spatial and Temporal Land-Use Changes and Scenario Simulations

In this study, the PLUS model was used to simulate the land use changes for the year 2030 under different scenarios of the urban agglomeration occurring on the northern slopes of the Tianshan Mountains. The model contained two main modules based on the land-expansion-analysis strategy (LEAS) and a CA model based on multi-type stochastic patch seeding (CARS). The basic principle of the model is to extract part of the land-use expansion for each type of land-use between the two periods of land-use change, calculate the degree of influence of each driver on the expansion of different land-use types using the random forest algorithm, and obtain the development probability of each type of land-use. Based on this, the spatial distribution of land-use in the initial year was used as the simulation benchmark, and the development probability was used as the constraint to simulate the land-use-change pattern dynamically in time and space. The model combines the CA model with the patch-generation simulation strategy to improve the ability to simulate real landscape patterns, and therefore the model is widely used in the field of land-use simulation and prediction [28].

The process is as follows:

(1) Selecting the driver.

Taking into account the spatial and temporal accessibility and quantifiability of various types of data, a total of 12 natural conditions related to accessibility and socio-economic factors, such as elevation, slope, evaporation (ET), precipitation (Pre), population density, GDP, distance from highways, distance from provincial roads, distance from county roads, distance from township roads, distance from railways, and distance from river systems, were selected according to the scale characteristics and geographical features of the study area. Factors taken as predictor variables and policy and urban planning taken as scenario principles and constraints [29] were used for the simulations and predictions. At the same time, the LEAS module of the PLUS model was combined with the site-expansion factor-contribution analysis presented in Figure 3, which can effectively reflect the contribution of different influencing factors in each land category and, to a certain extent, the reasonableness of the selected influencing factors. In this study, the elevation and climate factors presented a greater influence on forested grassland and cultivated land than other factors, and the roads had the greatest influence on urban growth. Finally, future land-use simulations of the study area were performed based on the factors presented in Figure 4.

(2) Land-use area weighting-parameter setting.

Neighborhood weights can reflect the difficulty coefficient of the mutual conversion of each land type, which is mainly influenced by the reversibility of each land type change, with values between 0 and 1. The larger the value, the higher the stability and the lower the probability of conversion [30]. This study referred to and compared the results of previous studies on other similar areas, combined the proportion of the extended area and land-use change characteristics of each land-use type in the study area from 2010 to 2020, inferred the difficulty of mutual conversion through the proportion of the extended area, and debugged the model to obtain the neighborhood weight values of each land-use type in the study area (

Table 1. Field weight parameters of each land-use type in the urban agglomeration on the northern slopes of the Tianshan Mountains.

Land-Use Types	Cropland	Forest	Grassland	Waterbody	Built-Up Land	Unused Land
Neighborhood weight	0.8	0.5	0.8	0.5	0.8	1.0

).

(3) Future land-use simulation scenario setting.

This study referenced the Tianshan North Slope Urban Agglomeration Development Plan (2017–2030) and the previous relevant studies of the study area, and was based on the policies related to economic development and ecological protection in the study area, the requirements for sustainable development, and the various possibilities for its future development, based on the land-use change characteristics of the Tianshan north slope urban agglomeration in the past and its future regional spatial development plan. Three typical scenarios, natural development (S1), ecological protection (S2), and rapid economic development (S3), were set up from the perspectives of the urban expansion rate and ecological protection, respectively, to predict the number and spatial distribution of land-use types in the Tianshan north slope urban group in 2030. Under the natural development scenario, the land use situation in 2030 was predicted by the Markov model based on the transfer probability matrix from 2010 to 2020; under the ecological protection scenario, reference was made to the “natural forest protection project” and ”basic farmland protection” policies, and the land-use situation in 2030 was also predicted by the Markov model. The ecological protection scenario was obtained by referring to the “natural forest protection project” and ”basic farmland protection” policies, adjusting them according to the historical natural development pattern; the rapid economic scenario was obtained by referring to the urban planning and development projections of the Tianshan north slope urban agglomerations.

(4) Simulation accuracy validation.

Using the PLUS model, the 2020 land-use map of the study area was simulated based on the 2010 land-use model and compared with the actual land-use status in 2020. As shown in Figure 5, the 2020 land-use map based on the 2010 simulation differs to some extent from the actual land-use map in 2020 as changes in the policy of returning farmland to forest and grass had affected the regional landscape of the study area. To verify accuracy of these findings, Kappa indices were calculated for both periods using the validation module of the PLUS model. The average Kappa coefficient for the six land-use types was 0.864, with an overall accuracy of 0.918, indicating that the accuracy of the model simulation can meet the accuracy requirements.

3.2.2. Analysis of Land-Use Intensity Change

Land-use intensity is closely related to ecosystems and is an important factor contributing to changes in ecosystem services. The equations are as follows [31] (Equation (1)). In this study, the land-use intensity of different land types was classified with reference to previous studies on land-use intensity indices [32,33,34] (

Table 2. Intensity index of land-use types in urban agglomeration on the northern slope of Tianshan Mountains.

Land-Use Types	Cropland	Forest	Grassland	Waterbody	Built-Up Land	Unused Land
Power level	3.0	2.0			4.0	1.0

):

$$L=100\times \left({\displaystyle \sum _{i=1}^{n}Pi\times Qi}\right)$$

(1)

where L is the composite index of land-use intensity; Pi is the land-use intensity of class i; and Qi is the area share of class i.

3.2.3. Assessment of Water Production

The principle of the water supply module is based on a water balance approach that takes into account factors such as temperature, vegetation, soil, and topography. In Equation (2), it is the amount of water produced per unit area by a raster cell over a certain period of time, i.e., it is the difference between the amount of precipitation per unit area by the raster cell and the actual evapotranspiration that is subtracted to obtain the water supply for each raster cell (which includes the surface runoff, the flow in the soil, and the amount of water held by interception and deadfall generated by the canopy), which is modeled using the following equations [18] (Equations (2)–(6)):

$$Y_x=\left(1-\frac{AE{T}_x}{p_x}\right)\times p_x$$

(2)

where AETx is the actual evapotranspiration produced by each grid cell x (in mm); Yx is the annual water production produced on each grid cell x (in mm); Px is the average annual precipitation produced on each grid cell x (in mm).

$$\frac{AE{T}_{\left(x\right)}}{P_{\left(x\right)}}=1+\frac{PE{T}_{\left(x\right)}}{P_{\left(x\right)}}-{\left[1+{\left(\frac{PE{T}_{\left(x\right)}}{P_{\left(x\right)}}\right)}^w\right]}^{1/w}$$

(3)

$$PE{T}_{\left(x\right)}=K_c{}_{\left(x\right)}\ast E{T}_o{}_{\left(x\right)}$$

(4)

$$W_{\left(x\right)}=\frac{AW{C}_{\left(x\right)}\ast Z}{P_{\left(x\right)}}+1.25$$

(5)

$$AW{C}_{\left(x\right)}=\min (\max Soi{l}_x,RootDept{h}_x)\ast PAW{C}_x$$

(6)

where PET(x) represents the actual potential evapotranspiration produced by each raster cell x (in mm); Kc(x) represents the reference evapotranspiration coefficient for each crop type; ET0 (x) represents the evapotranspiration produced by each vegetation type (in mm); W(x) represents the empirical non-physical parameter; AWC(x) represents the available water content for each vegetation type; Z represents the Zhang coefficient; Root Depth(x) represents the effective root depth for each vegetation type; Max Soil Depth(x) represents the maximum soil profile depth for each vegetation type.

The model sets parameters for different root depths based on different land-use data and the “Guide to calculating crop evapotranspiration—crop water requirements”; reference values for evapotranspiration coefficients (Kc) were set according to the Food and Agriculture Organization of the United Nations (FAO) to obtain the biophysical parameters needed for this study. The water yield assessment and validation process is as follows: by comparing with the water resources bulletin, repeated experiments are conducted to obtain the simulated value closest to the real value, and then the final input parameters of the model are obtained based on the results. The government’s water resources bulletin shows that in 2010 the surface water resources in Xinjiang were 1063 × 10⁶ m³. The modal water yield was calculated from the water resources bulletin, and when Z = 1, the modeled value for 2010 was 891 × 10⁶ m³, which is the closest to the actual value. Therefore, this study uses this value as the basis for the next step of the analysis.

3.2.4. Coordination Degree Model

To reflect the level of coordination between land-use intensity and water production, this study applied a coordination model to analyze the degree of coordination between the two variables, with the coordination degree taking on a value range of 0 to 1; the closer the value is to 1, the better the degree of coupling coordination. Firstly, land-use intensity and water production data were normalized and gridded, the study area was divided into 1 km × 1 km grids, and the mean values of land-use intensity and water production in each grid cell were extracted to calculate the coordination of the average annual growth rate between the two factors (

Table 3. Classification of coordination types.

Coordination Range	Category Division	Coordination-Type Area	Description
[0, 0.5)	Severely dysfunctional	Dysregulated decline	Land-use intensity and water production growth in a state of dysfunctional decline
[0.5, 0.6)	Mildly disordered
[0.6, 0.7)	face dissonance	Run-in coordination	Land-use intensity and water production growth in harmony
[0.7, 0.8)	Barely coordinated
[0.8, 0.9)	Basic coordination	Coordinated development	Represents a coordinated and rapid increase in land-use intensity and water production
[0.9, 1)	Good coordination

) [35]. The calculation was done as follows [36] (Equation (7)):

$$O=\frac{\left|\left(x+y\right)÷\sqrt{2}\right|}{\sqrt{x^2+y^2}}$$

(7)

where O is the coefficient of the coordination of changes in land-use intensity and water production; x and y are the average annual growth rates of land-use intensity and water production, respectively. The breakdown of the data used in this study is shown in

Table 4. Data sources and notes for this study.

No.	Datasets	Data	Data Resources
1	Land-use/cover datasets	Land-use/cover data	http://www.resdc.cn/data.aspx?DATAID=252 (accessed on 20 January 2021)
2	Vector datasets	Administrative district boundaries	https://xinjiang.tianditu.gov.cn (accessed on 15 March 2022)
3	Meteorological data sets	Average annual precipitation	http://data.tpdc.ac.cn (accessed on 15 March 2022)
		Average annual temperature
4	Social datasets	Population	http://www.resdc.cn/data.aspx?DATAID=252 (accessed on 20 January 2021)
		Gross domestic product	http://www.resdc.cn/data.aspx?DATAID=252 (accessed on 20 January 2021)
5	Road network data sets	Distance from railways	https://www.webmap.cn/main.do?method=index (accessed on 15 March 2022) calculated by Euclidean distance
		Distance from motorways
		Distance from provincial roads
		Distance from county roads
		Distance to country roads
		Distance to water
6	Topographic factor dataset	Elevation	http://www.gscloud.cn/#page1/2- (accessed on 20 January 2021)
		Slope
		Slope direction
7	Impact factor dataset	Nighttime lighting index	GEE (accessed on 15 March 2022)
		Normalized difference vegetation index	MODISTerra vegetation index data MYD13A1 16-day products
		Normalized difference building soil index	MODISA qua surface reflectance data MYD09A1 8-day products
8	Soil datasets	Clay, silt, sand	HWSD v1.2

.

4. Results and Analysis

4.1. Analysis of Spatial and Temporal Changes in Land-Use and Projections for the Study Area

The accuracy of the future land-use projection simulations has been verified to meet the accuracy requirements. The verification process is described in the methods. The land-use type of the urban agglomeration from 2000 to 2020 was dominated by unused land, which accounted for more than 50% of the total area, followed by arable land and grassland. The area of water bodies, forest land, and construction land only accounted for a small part of the study area (Figure 6). The total area of woodland and water bodies had decreased, with the area of woodland decreasing by 2526 km² and the share of water bodies decreasing by approximately 47.52% over a 20-year period. The grassland showed a decreasing trend followed by an increase, probably due to the fact that during the period 2000 to 2010, a large area of grassland in the central plains of the study area was converted for the purpose of economic development during the implementation of the “Western Development” strategy and the “One Belt, One Road” initiative. This had a significant impact on the local economy and simultaneously had a particular impact on the environment, which was later protected and restored through the relevant policies of “returning forests to grass and grazing” and environmental protection regulations. The largest increase in the area of land used for construction was from 0.83% to 1.77% of the original area (

Table 5. Land-use area from 2000 to 2020 and predicted area under three different in 2030 (km²).

	2000	2010	2020	2030
				S1	S2	S3
Land-use types	15,387	21,386	22,542	22,183	22,345	22,541
Forest	5139	2763	2613	2607	2604	2609
Grassland	62,325	58,024	59,072	59,676	58,934	59,676
Water body	3338	1851	1762	1639	1668	1767
Built-up land	1611	2346	3432	4439	3533	4439
Unused land	106,428	107,861	104,825	102,427	103,887	101,939

).

In the S1 scenario, the area of grassland and built-up land will increase by 604 km² and 1007 km², respectively, in 2030, compared to 2020, while the area of other land-use types will decrease. The area of arable land, water, and unused forest will decrease by 359 km², 123 km², and 2398 km², respectively, while the area of forest land will hardly change. The S2 scenario is based on the historical land-use change pattern of urban agglomeration and ecological protection as the dominant factor; therefore, land-use in this scenario will not change much compared to 2020. In the S3 scenario, the areas of forest land, grassland, and water will be 2609 km², 59,676 km² and 1767 km², respectively, and the area of unused land will change the most, from 104,825 km² in 2020 to 101,939 km² in 2030, while other land-use areas will change less.

4.2. Analysis of Land-Use-Intensity Changes in the Study Area

The most direct manifestation of the difference in the degree of human input and the use of different land resources is the difference in land-use intensity. In this study, according to the size of the study area and the calculation volume, the comprehensive land-use-intensity index was calculated separately for each grid, using 10 km × 10 km as a cell. The reasons for this spatial distribution are the following: the relatively high population density in these areas, the concentrated distribution of urban residential and arable lands, and the high interference of human activities, as well as the high concentration of human settlements and industrial activities. The intensity of land-use is relatively high in areas with the highest concentration of human and industrial activities. In terms of the general trends (Figure 7), land-use intensity from 2000 to 2030 is accompanied by an expansion of urban and rural construction land and an expansion of the high-intensity zone, with the area of the high-intensity zone being smaller in the S2 scenario than in the other two scenarios.

4.3. Analysis of the Change in Water Production of the Study Area

As can be observed in Figure 8, the water production in the Tianshan north slope urban agglomeration in 2000, 2010, and 2020 is 517.26 × 10⁶ m³, 582.28 × 10⁶ m³, and 456.38 × 10⁶ m³, respectively, showing an increasing trend followed by a decreasing trend; water production decreased by a total of 60.88 × 10⁶ m³ in 20 years. The water production values in scenarios S1, S2, and S3 in 2030 are 509.10 × 10⁶ m³, 510.90 × 10⁶ m³³, and 502.28 × 10⁶ m³, respectively. Compared to 2020, water production in 2030 will increase by 52.72 × 10⁶ m³, 54.52 × 10⁶ m³, and 45.90 × 10⁶ m³ in the three scenarios, respectively. Spatially, this presents a trend of high water production rates in the east, low water production rates in the west, high water production rates occurring again in the north, and low water production rates in the south, with the water-production function decreasing from the foothills of the Tianshan Mountains to the north and south during the study time period. In terms of land-use types, the contribution of different land-use types to the total regional water production rate is in descending order: grassland, cropland, forest land, construction land, and unused land. The area of grassland accounts for approximately 32% of the total area of the region and produces approximately 40% of the total water in the region, with a decreasing trend in the area but still increasing in water production, making it the most important water producing area in the region. Due to the small distribution of water in the study area, the area was one of high water production rates during the study time period. These trends continue over time from 2020 to 2030, with a decreasing trend in water production in the arable, forest, and unused lands and an increasing trend in water production in the grassland area in the three scenarios for 2020–2030. In the S2 scenario, the total water yield in the study area increased by 54.52 × 10⁶ m³.

The water production of cities within the Tianshan North Slope urban agglomeration also varies between 2000 and 2030, as shown in Figure 9. The high values of water production are mainly found in the western and north-central regions, with Shihezi and Fukang having higher water production in 2000 and Kuitun, Huyanghe, and Dushanzi having more than 32% of the water production by 2010. The lower values of water production were mainly in the northern and eastern areas of Shanshan and Tokson, where the average annual water production was around 16%. This distribution is related to the special geographical location of the study area.

4.4. Analysis of Changes in the Coordination between Land-Use Intensity and Water Production

From an overall perspective (Figure 10), the coordination type accounted for the largest proportion of the three scenarios, with an increasing then decreasing trend, accounting for 60%, 89%, and 87% in the periods 2000–2010, 2010–2020, and 2000–2020, respectively. From 2000 to 2010, the area of high land-use intensity in the study area expanded, land-use intensity was rapidly increasing at a quicker rate than water production, and both land-use intensity and water production rates significantly increased. From 2010 to 2020, land-use intensity stabilized, and the main source land types of regional water production changed relatively slightly, presenting a trend of coordinated over-expansion. The trends of the various types of coordination for the three scenarios in 2030 were roughly the same as those in the previous time periods, and the proportion of coordination occurring in the S2 scenario was better compared to that of the other scenarios. Due to the complex topography of the study area, the distance of forest land, the water areas, and unused land from built-up areas, the difficulty of development and construction, and the low level of human activities, the changes in land-use intensity and water production in the distribution areas of these land-use types were relatively limited; thus, the growth rates of both were basically synergistic and the coordination types were not prone to change, while unused land, as the main land-use type in the study area, accounted for approximately 55% of the study area. The coordination between land-use intensity and water production rates therefore present a more stable characteristic overall over the study time period.

4.5. Analysis of Factors Influencing Changes in Water Production Rates in the Study Area

This study was based on the spatial and temporal distributions of water production in the Tianshan north slope urban agglomeration and its variation patterns, combining natural and anthropogenic factors for quantitative research [37]. Using the data for 2020 as an example, the natural and anthropogenic data were gridded using the grid method to create 2098 grid points at a size of 10 km × 10 km within the study area, and were used to match the water yield (detection factor) with influencing factors (explanation factor). The interaction between the various factors and water production rates in the urban agglomeration of the northern slopes of the Tianshan Mountains, and whether there were significant differences between them, thus revealed the main driving forces influencing the spatial distribution.

As can be observed in Figure 10, the different factors are ranked in descending order of q-values in terms of their determinants of water production: NLI (0.38) > LUCC (0.23) > DEM (0.17) > AAE (0.16) > APP (0.14) > slope (0.11) > GDP (0.09) > NDVI (0.08) > NDBSI (0.06) > POP (0.02). From the results of the factor detection, we can observe that the first dominant factor affecting the spatial variation of water production is NLI, with a contribution rate of 0.38, followed by LUCC, DEM, AAE, and other factors with equal influences, while the remaining factors present a relatively small contribution rate to the spatial variation of water production.

In terms of the magnitude of the interaction (Figure 11), this can be observed by the following values: DEM ∩ AAP (0.33), AAP ∩ GDP (0.30), AAE ∩ AAP (0.28), DEM ∩ GDP (0.27), DEM ∩ LUCC (0.26), slope ∩ LUCC (0.24), AAP ∩ LUCC (0.21), and AAE ∩ LUCC (0.20). The interaction between the other factors and land-use types in the study area was significant, and the interaction was more pronounced than the interaction between other factors as different land-use type structures determine the distribution patterns of different ecosystem types, making the interaction between the LUCC factor and other factors more pronounced than other factors’ interactions, with the distribution of land-use types playing a major role in the changes in water production rates. The interaction of DEM, AAP, and AAE with other factors was also obvious, indicating that the topography and climatic conditions of the study area had a certain influence on the water yield in the study area.

5. Discussion

5.1. Future Projections of Land-Use Changes

Land is the basis for human survival and social development, and rational land-use is necessary to achieve the goal of sustainable urban development [38]. The ecological changes and their characteristics resulting from the interaction between human activities and the natural environment can be reflected more clearly through land-use changes [39]. In this study, from 2000 to 2020, arable and construction lands in the study area presented an increasing trend, forest land and water areas presented a decreasing trend, and grassland and unused land presented a decreasing trend followed by an increasing trend. In the early years of the “Western Development” strategy, a large number of new industrial and commercial areas, excessive groundwater exploitation, and overgrazing led to the natural degradation of forest and grass vegetation and a reduction in land coverage. This is consistent with the results obtained by the majority of scholars [40]. For the present study, area land-use changes occurring in 2030 under different scenarios were mainly dominated by construction-based land expansion, arable land degradation and grassland expansion, but the degree of change varied depending on the scenario. The economic development scenario presented the largest increase in built-up land area, while the ecological conservation scenario exhibited the largest increase in forest land area. However, it was predicted that such results may depend on the impact factor, and we therefore recommend that similar studies are conducted by researchers.

In simulating and predicting future land-use patterns, this study selected the PLUS model, which applies a new analysis strategy to better explore the causal factors of various types of land-use change. This strategy combines the advantages of the existing TAS and PAS, avoiding the need to analyze the types of transformations that exponentially increase alongside the number of categories and retaining the model’s ability to analyze the mechanisms of land-use change over a certain period of time. It helped us to achieve high-accuracy simulation data in the shortest period of time [28]. The PLUS model parameters’ transition matrix and neighborhood weights were used in the calibration process. The transtion matrix was adjusted empirically. Neighborhood weights were calculated as a percentage of the land-use extension for each land-use type. The above parameters were used to improve the accuracy of the simulations during the model calibration process. The increased intensity of human activity led to increased ecological pressures with serious implications for the global ecosystem [41,42,43]. Accurately predicting future urban development patterns is vital for managers and policy makers. The scientific validity of the data is linked to the techniques used. This study reduced the objectivity of the modelling results by considering the contribution of influencing factors.

5.2. Variation in Water Production and Its Influencing Factors

From 2000 to 2030, the water production rate in the urban agglomeration on the northern slopes of the Tianshan Mountains presents an annual decrease and a trend of high values in the east, low values in the west, high values in the north, and low values in the south, with the water production function decreasing from the foothills of the Tianshan Mountains to the north and south during the study period. Forests and grasslands are the most important water-producing areas in the region and play a more important role than other lands do in the study area. According to the analysis of the impact of land-use change on water production rates, it can be observed from the results that the change in water production rates was mainly due to the change in grassland and water areas, and that the change in land-use intensity and water production rates from 2000 to 2030 are in harmony. This conclusion can be verified further by the results of the analysis of the coordination between land-use change and water production rates from 2000 to 2030. The results of this study, combined with the analysis of land-use development and water yield changes in the study area, suggest that land-use change can have a positive impact on the ecosystem’s water yield to a certain extent, and that for the study area, the ecological conservation development scenario may be more in line with the future urban development model, with the construction of a good ecological environment as a basic guarantee for the development of urban agglomerations, which is conducive to the comprehensive and sustainable development of urban agglomerations. Our approach helps to identify hotspots of gains and losses in the ecosystem of water yield functions by using two simple models to assess the relative importance and combined impacts of land-use change on water yield functions, which can be used for more informed environmental investment decisions. Based on our results, a pronounced impact of land-use change on the urban scale can be identified. Our results suggest that wise land-use management can improve ecosystem services.

5.3. Research Limitations and Future Studies

Predicting the changes in an ecosystem’s water yield provides a novel way of observing the evolution of land-use change and of water yield in the future. As described in the methodological framework (Figure 2), the main elements of our approach were factor selection, water production function assessment, and factor analysis. Our method provides a straightforward way to explore the possible impacts of land-use change on ecosystem services. The method is simple and can be used in other regions. In this study, the impact of the changes in land-use intensity on water yield was taken into account to further validate the results of the study and to produce more reasonable and scientific conclusions. The PLUS and InVEST models could clearly reflect the future changes in the water-yield values in the region, and the prediction of the changes in the ecosystem’s water yield figures provided a novel way of studying the evolution of land-use and of water yield in the future. However, several uncertainties still remained in the study. Firstly, the method used to calculate future land demand in a more reasonable way by integrating various natural and social factors is one of the key ways to improving the accuracy of future land-use pattern predictions. Secondly, the InVEST model’s results are highly sensitive to the input parameters. As the Z parameter is the main input in the model, it must be set in conjunction with the measured data. In addition to the fact that land-use change is the dominant factor that influences water production, the present study verified this conclusion through the use of a coordinated model, and further research is required to achieve an insight into the model’s performance following land-use change. The use of a large amount of measured data to obtain dynamic water-yield coefficients can improve the accuracy of assessments of an ecosystem’s water yield, and the continuous monitoring of sample sites in the study area should be performed for many years to verify the reasonableness of the water yield coefficients.

6. Conclusions

This study is based on land-use modeling techniques and a coordination degree model, which not only helps to provide insight into future urban growth patterns in the Tianshan North Slope urban agglomeration but also helps to inform different land-use stakeholders. The study provides a spatial assessment of the water production function of ecosystems in arid and semi-arid regions, which helps in understanding the impact of land-use change on the water production function of the urban agglomeration on the northern slopes of the Tianshan Mountains. The combined impacts and relative importance of land-use change in terms of changes in ecosystem services were identified. A feasible and replicable framework was proposed for land-use decision-making in ecologically fragile areas. The results show that, at the scale of urban agglomerations, land-use change has a greater impact on water-producing functions than climate change does. In most of the areas we studied, land-use change had a dampening effect on water-producing functions at the urban scale. Our study also supports the hypothesis that rational land-use will, to some extent, have a positive impact on the water-producing function of ecosystems and that effective management can help managers to develop more comprehensive spatial plans. A combination of PLUS model scenario analysis, geo-probing, and consistency analysis can assess the impact of drivers on ecosystem services. Despite the limitations of this study, the findings are practical, methodological, and policy-relevant, and can support the use of ecosystem service information in land planning and the development of more effective ecosystem conservation strategies. Finally, this study in can serve as a reference for sustainable ecosystem services planning in the Tianshan North Slope urban agglomeration and other similar areas.