Evaluating Trade-Off and Synergies of Ecosystem Services Values of a Representative Resources-Based Urban Ecosystem: A Coupled Modeling Framework Applied to Panzhihua City, China

: Following signiﬁcant urban expansion, the ecological problems of resource-based cities are gradually exposed. It is of great signiﬁcance to study the ecosystem services of resource-based cities to achieve their sustainable development goals and to alleviate the conﬂicts between environmental protection and the utilization of the surrounding resources. However, in the current research on resource-based cities, few scholars have combined multiple minerals and multiple ecosystem services to explore the impact of mineral resources on the ecosystem. In this study, based on the historical data spanning from 2002 to 2018, we used the CA–Markov model to project the land use of Panzhihua City to 2030. Based on future land use projection, we quantiﬁed four ecosystem services (ESs) variables, including water yield, carbon storage, habitat quality, and soil conservation, using the InVEST model from the perspective of land use evolution in Panzhihua City. In addition, we explored the trade-offs and synergies of different ecosystem services and the correlations between different mineral species and ecosystem services using Spearman’s correlation coefﬁcient. Results showed the following: (1) During 2002–2018, water yield service, habitat quality service, and carbon storage service of Panzhihua City decreased year by year, and soil conservation service showed signiﬁcant ﬂuctuations; most of the low ESs areas were distributed in the central region of Panzhihua. On the contrary, most high ESs areas were located in the forest region. (2) The trade-offs and synergistic relationships among different ecosystem services showed signiﬁcant spatial variations. There were synergistic relationships among ESs and weak trade-offs between water yield services, soil conservation, and habitat quality services. There was also signiﬁcant spatial variability in the trade-offs and synergies among ecosystem services, with water production services showing “east trade-offs and west synergies” with soil conservation and habitat quality services, and most of the rest showing trade-offs in urban areas. (3) ESs in mining areas showed trade-offs in general, mainly between water production services and carbon storage services, with clay as the major negative factor of mineral species, and iron ore mines that have undergone ecological protection construction showed the lowest negative impact on ecology.


Introduction
Ecosystem service value assessment is an essential reference and basis for ecological and environmental protection, ecological function zoning, environmental and economic accounting, and ecological compensation decision-making [1,2]. With the acceleration of economic development and urbanization, the decline of ecosystem services and the decay of biodiversity have become serious global problems. The assessment of ecosystem service functions helps provide a basis for the rational use of resources and then contributes to the sustainable development of the ecosystem.
Since the end of the 20th century, with the gradual clarity of the concept of ecosystem services and the deepening of the research on the division of ecosystem services, the research on ecosystem services has made great progress [3,4]. Domestic scholars have classified the indirect value of ecosystem service functions from the ecology perspective and constructed a value equivalent factor table based on the previous work, which provides a theoretical basis and reference for natural property assessment and ecological compensation [5,6]. With the deepening of research, evaluating ecosystem service value has become a research hotspot. Scholars have researched the methods of physical quality evaluation and value quantity evaluation of ecosystem services [7,8].
The assessment of physical quality methods of ESs is usually based on ecosystem service processes or ecosystem service functions, which can objectively reflect the formation mechanism of ecosystem services [9][10][11]. Strengthening the integration of ecosystem processes and services is the current trend in ecosystem assessment. With the deepening of ecosystem services research, various assessment models have emerged. Currently, the more common and open source models used internationally are InVEST model, the ARIES model, the SoLVES model, etc. Among them, the ARIES model is currently only applicable to some case areas in the United States and cannot take ecological or socioeconomic influences into account when used in other areas, but it has high accuracy for assessing functional areas and has excellent potential for development. The SoLVES model assesses ecosystem service functions innovatively, such as aesthetics, biodiversity, and recreational life in terms of public attitudes and preferences. However, when it is applied to agroecosystems, there is the problem of using the same landscape parameters for different landscape types. Although some data are difficult to collect, the InVEST model has strong applicability to small amounts of demand data. Furthermore, all modules of it are independent, allowing users to input relevant data of the study area, so it is suitable for a wide range of research [12,13]. Currently, the InVEST model has more applications for water yield, biodiversity, carbon storage, and soil erosion modules and fewer applications related to pollination and aesthetic assessment modules. Both domestic and international studies tend to focus on individual ecosystem services [14][15][16], the impact of land use change on ecosystem services [17][18][19][20], and the analysis of spatial and temporal patterns of ecosystem services [21,22], among which water yield services have been studied with more enthusiasm.
Land use change is a complex dynamic system with the characteristics of change discontinuity, landscape matrix, mixed land use categories, irreversible change, etc. [23][24][25][26]. Among the relevant models for land use simulation and prediction, it is challenging to predict the spatial pattern changes of land use with the traditional Markov model [27]. Many scholars internationally have conducted studies on urban growth using the Cellular Automata (CA) model [28][29][30][31], which is highly capable of simulating the spatial and temporal evolution of spatially complex systems, but due to the local interactions between system elements and the single control element, it is difficult to reflect the social, economic, and other macro factors affecting the regional ecological security pattern with this model. Some scholars used the CA-Markov model to simulate and predict land use changes, and the studies show that the current land use patterns still have many urgent problems to be regulated [32,33]. The CA-Markov model integrates the ability of the CA model to simulate the spatial change of complex systems with the advantages of the Markov model in long-term prediction. It not only improves the prediction accuracy of land use type transformations but also effectively simulates the spatial change of land use patterns, which is scientific and practical [34,35].
When people are involved in land use management, they often pursue only one or several types of ESs [36]. However, this method intentionally or unintentionally impacts the provision of other ecosystem services, thus raising the issue of trade-offs and synergies of ecosystem services. Ecosystem service trade-offs are situations in which the use of one ecosystem service increases or causes another benefit to decrease; synergies are situations in which two ecosystem services increase or decrease simultaneously [37]. Neglecting ecosystem service trade-offs or synergistic relationships may reduce the supply capacity of specific ecosystem services [13]. Among the investigations of ecosystem service trade-offs and synergies, scholars tend to focus on regional ecosystem services [38][39][40][41], the influence of natural resources on ecosystem service trade-offs and synergies [42][43][44], the influence of land use change on ecosystem services [44,45], and the influence of urban spatial planning on ecosystem services [46], with fewer studies combining mineral resource development and resource-based urban ecosystem service trade-offs and synergies. In the research on ecosystem services of resource-based cities, domestic scholars constructed a composite urban carrying capacity index based on mineral resources from the synergy of economy, society, resources, and environment, which provides new perspectives for investigating the carrying capacity of mining-based cities [47]. Generally, the research content is biased towards the impact of land use change on ecosystem services and ecosystem vulnerability [48][49][50][51]. Few related studies explore the synergistic trade-off relationship between ecosystem services and mineral species. The irrational use of resources in resource-based urban ecosystems leads to excessive environmental depletion, and the declining function of natural ecosystems is becoming increasingly severe [52,53]. Therefore, it is desirable to investigate the ecosystem service functions of resource-based cities to help explore the sustainable development path of resource-based cities [54]. In the research of combined mineral resources, some scholars have explored the impact of mineral resources on the ecosystem from the perspective of single mineral species such as a coal mine [55,56] or a single ecosystem service [57]. Multiple ecosystem services and multiple minerals are rarely considered. In this study, we investigated the impacts of multiple mineral species on ecosystem services and explored the trade-offs and synergy between different mineral species and ecosystem services.
Panzhihua is a pivotal city for mineral resource development in China, so it is an important topic to investigate the relationship between different mineral species and ecosystem service functions in Panzhihua. The adjustment plan of Panzhihua's overall land use plan proposes to strengthen the land use planning and management and to transform the city from an industrial and mining base to an ecologically livable city. The "14th Five-Year Plan for Ecological and environmental protection of Panzhihua" pointed out that it is necessary to accelerate the development of low-carbon industries and to systematically promote ecological and environmentally friendly construction. In this paper, based on the land use data of Panzhihua city in 2002,2006,2010,2014, and 2018, we used the CA-Markov and InVEST models to assess and predict the evolutionary characteristics of the ecosystem services in the Panzhihua area from the perspective of land use. Furthermore, we used the Spearman correlation coefficient to analyze the trade-off and synergistic effects between different mineral species and ecosystem services and thoroughly combined it with the distribution of mineral resources in Panzhihua city to explore the correlation between mineral species and ecosystem services to provide a reference for the sustainable development and rational planning of resource-based cities.

Study Area
Panzhihua (26 • 05 -27 • 21 N, 101 • 08 -102 • 15 E) is located in the southernmost part of Sichuan Province, China, in the middle and southern part of the Panxi Rift Valley, with many hills and mountain plain canyon landscapes, and the terrain sloping from northwest to southeast (Figure 1). The climate is a compound of various climate types, such as south subtropical and north temperate, with an average annual temperature of about 20.8 • C, yearly precipitation of about 950 mm, distinct dry and rainy seasons, and a significant difference in temperatures between day and night. The city's total area is 7440.398 km 2 , with a resident population of more than 1.2 million. Panzhihua, as one of the four major iron ore regions in China, is rich in mineral resources. Seventy-six types of minerals have been discovered, and 7.18 billion tons of iron ore (mainly vanadium and titanium magnetite) have been proven. The associated vanadium, titanium, graphite, and other resources are at the forefront in China. Cobalt, chromium and other rare metal minerals, and coal, dolomite, granite, and other non-metallic minerals are abundant. development and rational planning of resource-based cities.

Study Area
Panzhihua (26°05′-27°21′N, 101°08′-102°15′E) is located in the southernmost part of Sichuan Province, China, in the middle and southern part of the Panxi Rift Valley, with many hills and mountain plain canyon landscapes, and the terrain sloping from northwest to southeast (Figure 1). The climate is a compound of various climate types, such as south subtropical and north temperate, with an average annual temperature of about 20.8 °C, yearly precipitation of about 950 mm, distinct dry and rainy seasons, and a significant difference in temperatures between day and night. The city's total area is 7440.398 km 2 , with a resident population of more than 1.2 million. Panzhihua, as one of the four major iron ore regions in China, is rich in mineral resources. Seventy-six types of minerals have been discovered, and 7.18 billion tons of iron ore (mainly vanadium and titanium magnetite) have been proven. The associated vanadium, titanium, graphite, and other resources are at the forefront in China. Cobalt, chromium and other rare metal minerals, and coal, dolomite, granite, and other non-metallic minerals are abundant.

Data Sources
Multiple sources of data were used in this study. For more information on the data and data sources, see Table 1. We preprocessed the data in PIE-Basic remote sensing image processing software. We reclassified the land use data into seven land classes: crop, forest, shrub, grassland, water, barren, and impervious, and we coded these classes from 1 to 7, respectively. We extracted slope rasters based on DEM data and extracted rivers and watersheds through hydrological analysis tools in ArcGIS. We used ArcGIS to calculate Euclidean distances to obtain grids of distances to rivers, roads, railways, and settlements. We screened 20 meteorological stations in Panzhihua City and surrounding areas and obtained daily rainfall and annual average rainfall grids by inverse distance weighted interpolation in ArcGIS. The inverse distance weighted method takes the distance between the interpolation point and the sample point as the weighted average weight. The sample points closer to the interpolation point are given more weight. All spatial data were unified under the Asia_North_Albers_Equal_Area_Conic projection coordinate system, and the raster data were resampled to 30 m spatial resolution.

Methodology
We innovatively combined the CA-Markov Forecasting model with the InVEST Forecasting model in the study. We forecasted land use in 2030 and selected four modules, namely, water yield, soil conservation, habitat quality, and carbon storage, to evaluate the value of ecosystem services and to construct an ecosystem services evaluation system for resource-based cities. In particular, we optimized and improved the evaluation unit of the mining ecosystem. We constructed a new model of ecosystem service trade-offs and synergies across different mining areas according to the mining type ( Figure 2). The cellular automata (CA) model is a dynamic model that presents discrete data in time, space, and state, constructed from temporal causality and spatial interactions, and it The cellular automata (CA) model is a dynamic model that presents discrete data in time, space, and state, constructed from temporal causality and spatial interactions, and it is able to model the spatio-temporal evolution of complex systems [60]. The following equation can express it: where S is the metacell space, which is the set of finite discrete metacells; f is the transition rule function of metacell states; N is the neighborhood of each metacell; and t + 1 and t are two different moments.
(2) Markov Model Markov is a special motion process based on the theory of the stochastic process of the mathematician Markov, which has "no after effect" and can collect complex information in the form of a state transfer matrix using historical factors and present states for statistical and developmental law exploration [61]. Since land use changes have similar properties, Markov models are widely used to predict dynamic changes in land use [62]. The Markov model is calculated using this equation.
where S (t + 1) and S (t) are the land use states at the moment t + 1 and the moment t, respectively, and P ij is the land use transfer probability matrix. (

3) CA-Markov Model
The CA-Markov model is a hybrid model consisting of cellular automata (CA) and a Markov chain. In this study, we used IDRISI Selva v.17 to predict the future land use of the study area on the CA-Markov model. IDRISI Selva v.17 is a platform that combines image processing and geographic information analysis [63,64]. We used the CA-Markov module in it to obtain the transition probability and the transition area matrix [65]. Land uses of 2014 and 2018 were used to construct the transition probability matrix for the purpose of producing a simulated land use map for 2030.

Ecosystem Services Assessment
(1) Water Yield The water yield module in the InVEST mode is based on the Budyko principle of coupled hydrothermal equilibrium [66], which combines the effects of factors such as spatial differences in soil permeability and evapotranspiration of different land use types on runoff to construct a suitable model and takes raster as a unit to quantitatively estimate water supply capacity [67].
where Y xi is the annual water yield on the raster x when the land use is type i (mm); P x is the average annual precipitation of the raster x (mm); and AET xi is the actual average annual evapotranspiration (mm). where PET xi is the annual average potential evapotranspiration on the raster x when the land use type is i (mm); w x is a non-physical parameter of natural climatic soil properties; AWC x is the effective soil water content of the raster x (mm); and Z is the Zhang coefficient, which is a seasonal constant taking values between 1 and 30. In this study, the data from hydrological stations in the study area and relevant papers were analyzed [68], and the Zhang coefficient was continuously adjusted to debug the optimal assessment results.
(2) Soil Conservation In the InVEST model, soil conservation (SC) is obtained by subtracting the actual soil erosion (USLE) under manual management and conservation measures from the potential soil erosion (RKLS) under natural vegetation protection.
where R is the rain fall erosivity, calculated using a monthly calculation formula [69]; K is the soil erodibility factor, calculated using the EPIC model and modified with the correction method [70,71]; LS is the slope length-gradient factor; C is a cover-management factor, which is the ratio of soil loss from vegetated land or managed fields to soil loss from recreational land with continuous light tillage under the same environmental conditions [72]; and P is the soil conservation measure factor, which is the ratio of soil loss from sloping land with soil conservation measures to soil loss from sloping land without any measures under the same environmental conditions [73,74]. (

3) Habitat Quality
This module of the InVEST model is mainly used to obtain the distribution of habitat quality and the distribution of degradation in the study area by the sensitivity of different land types to each threat source and habitat threat density data and to assess biodiversity by the level of habitat quality, which is calculated as follows: where Q xj denotes the habitat quality score of raster x in land use type j; H j denotes the habitat suitability; and K is the semi-saturation parameter. In this study, the InVEST model manual was used as the basis, and the threat factor and sensitivity were set according to Yubin Bao [19], Shu Feng [75], Dazhi Zhang [76], and Xiaoyu Niu [77].

(4) Carbon Storage and Sequestration
In the InVEST model, the carbon storage and sequestration of the ecosystem (C total ) consists of four basic carbon pools: Aboveground carbon stock (C above ), belowground carbon stock (C below ), soil carbon stock (C soil ), and dead organic carbon stock (C dead ), and the model is calculated as follows: Determining the carbon density value of each carbon pool is the key to calculating ecosystems' carbon stock functions. In this paper, based on the carbon density estimation results of Xie Xianli et al. and Li Ke-Jean et al. for different land use types in China [78,79], we corrected the carbon density by using the correction method proposed by Alam et al. in combination with the precipitation correction factor [80].

Correlation Analysis
Sampling points were evenly laid out in GIS software to extract the four ecosystem services during the study period . Then the Spearman correlation coefficient between the two ecosystem services was calculated using mathematical and statistical tools. Finally, the Spearman correlation coefficient was applied to calculate the correlation of non-normally distributed data, which is calculated as follows [81,82]: where r 12 is the correlation coefficient between the two ecosystem services; n is the number of sampling points; ES1 i is the rank of the i-th sampling point of the first ecosystem service; and ES2 i is the rank of the i-th sampling point of the second ecosystem service. The larger the absolute value of r 12 , the stronger the correlation [56]. If it is positive and passes the significance test, there is a synergistic relationship between the pair of ecosystem services. If it is negative and passes the significance test, there is a trade-off relationship between the pair of ecosystem services [83].

Spatial and Temporal Changes in Land Use
In this study, land use data for a total of five periods at each four-year interval from 2002 to 2018 were selected, the land use dataset was cropped according to the administrative boundaries of the study area, and a land use transfer matrix was constructed (Table A1). Based on this, a land use change chord diagram was introduced to visualize the land use transfer changes (Figure 3).
where r is the correlation coefficient between the two ecosystem services; n is the number of sampling points; ES1 is the rank of the i-th sampling point of the first ecosystem service; and ES2 is the rank of the i-th sampling point of the second ecosystem service. The larger the absolute value of r , the stronger the correlation [56]. If it is positive and passes the significance test, there is a synergistic relationship between the pair of ecosystem services. If it is negative and passes the significance test, there is a trade-off relationship between the pair of ecosystem services [83].

Spatial and Temporal Changes in Land Use
In this study, land use data for a total of five periods at each four-year interval from 2002 to 2018 were selected, the land use dataset was cropped according to the administrative boundaries of the study area, and a land use transfer matrix was constructed (Table  A1). Based on this, a land use change chord diagram was introduced to visualize the land use transfer changes ( Figure 3). From 2002 to 2018, forest was the dominant part of land use in Panzhihua, accounting for more than 60% of the total area; the second type was cropland, accounting for 26.69% of the total area in 2018; followed by grassland, accounting for 8.31% of the total area in 2002; the remaining land use types include impervious surfaces, shrubs, barren land, and water, which are less than 5%. Long-term observations of land use show that cropland, From 2002 to 2018, forest was the dominant part of land use in Panzhihua, accounting for more than 60% of the total area; the second type was cropland, accounting for 26.69% of the total area in 2018; followed by grassland, accounting for 8.31% of the total area in 2002; the remaining land use types include impervious surfaces, shrubs, barren land, and water, which are less than 5%. Long-term observations of land use show that cropland, water, and impervious surfaces in Panzhihua continue to increase, while barren, grassland, and forest continue to decrease. Influenced by urbanization and industrialization, impervious surfaces expanded significantly from 2006 to 2010, with an increase of 24.38%, mainly from the conversion of cropland and grassland. On the other hand, the strong advance in agricultural production has dramatically increased the cropland area and decreased forest and grassland, with up to 15,181 hm 2 converted from forest to cropland from 2010 to 2014. In addition, shrubs were mainly converted from forest and farmland, with less increase or loss of wasteland and watershed.

Land Use Forecast Analysis
In this study, the prediction of dynamic land use changes in 2030 was carried out with the help of IDRISI. A total of seven driving influence factors, namely, DEM, slope, distance from river network, distance from road, distance from the highway, distance from railroad, and distance from the settlement, was selected to construct the suitability atlas and combined with the transfer matrix to build 5 × 5 standard filters with four iterations for the 2018 land use prediction with a Kappa coefficient of 0.8802 [84,85], which meets the criteria for the 2030 prediction. In the 2030 projection results, there were significant growths in cropland, shrubs, grasslands, and impervious surfaces, among which shrubs grew at a rate of more than 0.75, and all other growth was around 0.16; there were significant decreases in forests and bare lands, among which forests decreased by about 600 km 2 , and bare lands decreased at a rate of more than −0.75; and waters remained stable with a slight increase.

Spatial and Temporal Analysis of Ecosystem Services
In this study, quantitatively, the InVEST model was used to evaluate four ecosystem services in the study area. The evaluation results of six periods from 2002 to 2030 were obtained, as shown in Figure A1. The annual water yield showed a general downward trend from 2002 to 2018 and will improve slightly by 2030. The water yield of the northeastern part of Panzhihua was the best, while the central and northwestern parts were poor. Soil conservation had a trend of increasing, then decreasing, then increasing again. The central, northern, and eastern edges of Panzhihua were along the high-value area of soil conservation services, while the southern and southeastern areas were low-value areas. Habitat quality in the study area decreased year by year. Habitat quality was better in the central and northern parts of Panzhihua, while it was worse in the south and northeast. Carbon sequestration showed a decreasing trend, with higher stocks in the mountainous areas in the north-central part and lower stocks in the south-central and southeastern parts with higher urbanization levels.
The overall pixel distribution, spatial distribution, and annual average change rate of ecosystem services in the study period showed that the long-term time series of the four ecosystem services in the study area changed significantly, and the spatial characteristics were very different (Figure 4). Figure 4a reveals a slight increase in median water production in 2030 compared to 2018, and Figure 4i shows a continued slowdown in the decline in water production from 2006 to 2018, which collectively indicates a gradual improvement in water production services, which is well reflected in the local magnification in Figure 4e. Figure 4b shows that soil retention averages are concentrated at lower levels for most rasters. There are localized differences in soil retention between 2018 and 2030 in terms of increases and decreases (Figure 4f). Soil retention mean values repeatedly fluctuated during 2002-2018, with average annual rates of change exceeding 10% several times (Figure 4j). Figure 4c displays that the habitat quality of most pixels is above 0.8, but there is an overall decline in 2030 compared to 2018, reflected in a decrease in the number of high habitat quality pixels and a general reduction in the median and interquartile range. Figure 4g presents a sharp drop in regional habitat quality, with the annual decline in habitat quality leading to the increase and then slowing down (Figure 4k). Figure 4d indicates decreasing carbon sequestration in high-value areas and increasing numbers in low-value areas in 2030, with a significant decrease in regional carbon sequestration (Figure 4h). The average rate of change in carbon sequestration was similar to that of habitat quality, which presented a continuous trend of increasing decline until 2018 (Figure 4l). habitat quality leading to the increase and then slowing down (Figure 4k). Figure 4d indicates decreasing carbon sequestration in high-value areas and increasing numbers in lowvalue areas in 2030, with a significant decrease in regional carbon sequestration ( Figure  4h). The average rate of change in carbon sequestration was similar to that of habitat quality, which presented a continuous trend of increasing decline until 2018 (Figure 4l).

Spatial and Temporal Analysis of Ecosystem Services
A 500m × 500m fishing net was set up, and the ecosystem services were processed using the linear function normalization method [86]. The central image element values of the net were extracted to calculate the correlation coefficients and obtain the correlation coefficient matrix ( Figure 5). The upper triangle of the matrix shows the correlation coefficients and significance between ecosystem services, the lower triangle shows the scatter plot of each pair of ecosystem services, and the main diagonal shows the density distribution curve.

Spatial and Temporal Analysis of Ecosystem Services
A 500m × 500m fishing net was set up, and the ecosystem services were processed using the linear function normalization method [86]. The central image element values of the net were extracted to calculate the correlation coefficients and obtain the correlation coefficient matrix ( Figure 5). The upper triangle of the matrix shows the correlation coefficients and significance between ecosystem services, the lower triangle shows the scatter plot of each pair of ecosystem services, and the main diagonal shows the density distribution curve. Remote Sens. 2022, 14, x FOR PEER REVIEW 13 of 27  HQ and CS were significantly positively correlated from 2002 to 2018 (0.600 < r < 1, p < 0.001), indicating that they have a significant synergistic relationship. Their R-values peaked in 2010 and then gradually decreased, and the R-value was predicted to drop to 0.587 in 2030 (p < 0.001). SC was positively correlated with HQ and CS, reflecting a synergistic relationship. However, the R-value of SC and HQ increased yearly after 2006, indicating that the synergistic relationship was strengthening while the synergistic relationship between SC and CS became weaker. WY and SC were weakly synergistic, and their R-values (0 < r < 0.200, p < 0.001) were significantly lower than those of the previous three. From 2002 to 2014, WY and HQ were negatively correlated, showing a weak trade-off relationship and the trade-off is predicted to disappear in 2030. WY and CS also showed a weak trade-off relationship. The results indicate that the ecosystem services in Panzhihua tend to be coordinated, and that Panzhihua has played a positive role in maintaining the stability of the ecosystem by promoting ecological and environmental protection and construction and the coordination of ecosystem services.

Spatial Analysis of Ecosystem Service Trade-Offs and Synergies
The study area was sampled raster-by-raster, correlation coefficients between ecosystem services were calculated at the image element scale, and trade-offs and synergistic relationships were mapped according to the following classes ( Figure 6): strong synergy (r > 0, p < 0.05); medium synergy (r > 0, 0.05 p < 0.1); weak synergy (r > 0, p ≥ 0.05); weak trade-off (r < 0, p ≥ 0.05); medium synergy (r < 0, 0.05 ≤ p< 0.1); and weak synergy (r < 0, p < 0.05). The results indicate that the ecosystem services in Panzhihua tend to be coordinated, and that Panzhihua has played a positive role in maintaining the stability of the ecosystem by promoting ecological and environmental protection and construction and the coordination of ecosystem services.

Spatial Analysis of Ecosystem Service Trade-Offs and Synergies
The study area was sampled raster-by-raster, correlation coefficients between ecosystem services were calculated at the image element scale, and trade-offs and synergistic relationships were mapped according to the following classes ( Figure 6): strong synergy (r > 0, p < 0.05); medium synergy (r > 0, 0.05 p < 0.1); weak synergy (r > 0, p ≥ 0.05); weak trade-off (r < 0, p ≥ 0.05); medium synergy (r < 0, 0.05 ≤ p<0.1); and weak synergy (r < 0, p < 0.05). The spatial distribution of the relationship between WY, HQ, and CS was similar, with the strong synergistic relationship concentrated in the central main urban area and the strong synergistic relationship scattered in the northwest, northeast, and southeast The spatial distribution of the relationship between WY, HQ, and CS was similar, with the strong synergistic relationship concentrated in the central main urban area and the strong synergistic relationship scattered in the northwest, northeast, and southeast mountainous areas. The weak synergistic relationship between SC and HQ appeared in the central and south-central areas and decreased in a radial pattern to the surrounding areas. A weak trade-off relationship appeared in the study area's northwestern, northeastern, and southern parts. There was large insignificance between SC and CS, and there was a fine-grained trade-off area in the central part of the study area.

Impact of Mineral Development on Ecosystem Services
Based on the vector mineral data, the spatial distribution of nine essential minerals in Panzhihua, including coal, iron ore, and saprolite, was determined by mapping synthesis (Figure 7). In addition, the correlation analysis method was extended to mining areas, and the correlation coefficients of ecosystem services in mining areas were calculated pixel by pixel according to mineral types (Figure 8).
Remote Sens. 2022, 14, x FOR PEER REVIEW 15 of 27 mountainous areas. The weak synergistic relationship between SC and HQ appeared in the central and south-central areas and decreased in a radial pattern to the surrounding areas. A weak trade-off relationship appeared in the study area's northwestern, northeastern, and southern parts. There was large insignificance between SC and CS, and there was a fine-grained trade-off area in the central part of the study area.

Impact of Mineral Development on Ecosystem Services
Based on the vector mineral data, the spatial distribution of nine essential minerals in Panzhihua, including coal, iron ore, and saprolite, was determined by mapping synthesis (Figure 7). In addition, the correlation analysis method was extended to mining areas, and the correlation coefficients of ecosystem services in mining areas were calculated pixel by pixel according to mineral types (Figure 8).  Overall, the mining ecosystem services were moving toward trade-offs. The study area had a slight trade-off between WY and CS. However, within the amphibolite, granite, and chert mines, the R-value between WY and CS was below −0.600 (p < 0.001), reflecting a significant trade-off. The R-values between WY and HQ were negative in six of the nine mineral zones, and their trade-offs were stronger than the whole study area. The relationship between SC and HQ was positively correlated within mineral zones, and its R-value was even higher than those of the entire study area in some mineral zones, indicating Overall, the mining ecosystem services were moving toward trade-offs. The study area had a slight trade-off between WY and CS. However, within the amphibolite, granite, and chert mines, the R-value between WY and CS was below −0.600 (p < 0.001), reflecting a significant trade-off. The R-values between WY and HQ were negative in six of the nine mineral zones, and their trade-offs were stronger than the whole study area. The relationship between SC and HQ was positively correlated within mineral zones, and its R-value was even higher than those of the entire study area in some mineral zones, indicating stronger synergistic relationships than those of the entire study area. HQ and CS were positively correlated within mineral zones, but their synergistic relationships were weaker than the whole study area. The relationship between WY and SC and that between SC and CS had different trade-off/synergistic relationships across the mine area and did not show consistency.
The varying impacts of different mineral extraction types on ecosystem services were reflected in the differences in the ecosystem service correlation coefficients of those in separate mining areas. Compared to other mining areas, the R-values between ecosystem services in coal and iron ore mines were the least different from those in the whole study area, indicating that they had the least impact on ecosystem services. The trade-off between WY and CS was the most significant in the granite and tuff mining areas. There was a strong trade-off relationship between SC and HQ in dolomite mines. In shale, clay, and construction sand ore zones, WY reflected a significant trade-off with SC, HQ, and CS, while SC reflected significant synergy with HQ, SC with CS, and HQ with CS. The clay mine area had the largest difference in R-value from the entire study area, indicating that its ecosystem services were most affected.

Discussion
This study assessed ecosystem services and trade-off and synergy analysis in the Panzhihua region and explored the impact of mineral extraction on ecosystem services, filling the research gap in ecological analysis and mineral impacts in resource-based cities.

Model Selection and Parameter Modification
In this paper, we assessed the ecosystem services in Panzhihua based on the InVEST model. The InVEST model has nearly 20 ecosystem service modules, including terrestrial, freshwater, and marine ecosystem service assessment models, covering many aspects of ecosystem services. Four of the most widely used and important modules were selected in this study, namely, water yield, soil conservation, habitat quality, and carbon storage, which is consistent with previous studies [87]. Considering the rapid population growth and the increasing demand for water by human activities, as well as the fact that Panzhihua has rich forest resources but serious soil erosion, we selected these four ecosystem services. Other ecosystem service assessment modules for the model will be carried out in the next step. We chose four ecosystem service assessment modules, namely, water yield, soil conservation, habitat quality, and carbon storage, to conduct the study because some of the modules in the InVEST model are not well applicable in the region. In comparison, these four modules have high research value due to their strong relevance to production and life, and the study by Xiuming Wang [88] and Conghong Huang et al. [89] proved the applicability and reliability of the combination of these four modules for the study, and the accuracy of ecosystem service assessment was initially improved through the screening of modules. Moreover, we modified the parameters of the InVEST model considering the regional characteristics of Panzhihua. The module of water yield was moderately corrected for the Zhang coefficient based on the statistical data of hydrological stations in Panzhihua, which is consistent with the correction method of Dou Miao [68] and Xu Jianning et al. [90]. The result was close to the observed average value of hydrological stations, where the Zhang coefficient was 3.2. The soil conservation module was made applicable to the region by the correction method proposed by Keli Zhang [71], which is consistent with that of Shuo Wang et al. [91]. The carbon storage module uses the precipitation correction factor combined with the method proposed by Alam et al. [80]. The correction coefficients for above-ground and below-ground carbon density were 1.7994, and for soil density they were 1.1979, which is consistent with the correction method of Xijin Ren et al. [92]. Through parameter modification, the research method was more suitable for the Panzhihua area, and the accuracy of ecosystem services assessment was further improved.

Assessment of Ecosystem Services and Exploration of Land Use Impacts
The results of the study indicated that the overall trend of the physical quality of ecosystem services in the Panzhihua region was decreasing, contrary to the effects of value quantity studies by Liu Sha [93] and Zhao Haifeng [94], but according to Zhao Jingzhu et al. [7] the use of two different methods, namely, physical quality and value quantity for service evaluation of the same ecosystem often leads to other or even opposite conclusions, providing support for the results of the study. Central Panzhihua, as a key region for economic development, formed a clear low-value zone for ecosystem services, which is consistent with the findings of Liu Ting [95] and Li Daoan [96], indicating that human activities such as mineral resource development and urbanization development affect the development of the ecological environment. Relevant measures are urgently needed to strengthen ecological restoration construction. However, the differences between the 2002 and 2030 ecosystem service assessment results were insignificant ( Figure A1), mainly due to the large study area, the absence of major land use type shifts during the study time frame, and the stable climate and vegetation conditions during this period [97]. Combining the spatial and temporal changes in land use in the Panzhihua area, it could be revealed that the overall decreasing trend of ecosystem service quality was related to the conversion of forests and grasslands to agricultural lands and impervious surfaces in the Panzhihua area (Figure 3), which is consistent with the study of Luo Jing et al. [98].
In addition, except for the weak trade-off relationship between water production services and soil conservation and habitat quality services, the ecosystem services in the Panzhihua area were mainly synergistic, which is consistent with the findings of Chen Xinmeng [99], Wang Xi [100], and others, which indicates that the ecosystem in Panzhihua area has better stability. For the weak trade-off relationship presented between water yield and soil conservation and habitat quality services, it is mainly due to the richness of regional vegetation such as forests, which have high carbon storage and soil conservation services, and the high evapotranspiration of vegetation, which leads to low water yield services. For this situation, the land use types are reasonably adjusted according to the different ecological needs of different areas, such as steep slope areas, considering appropriate discarding of water-producing ecosystem services, and ensuring soil conservation services, thus preventing excessive soil nutrient loss. The Panzhihua municipal government proposed promotion of the construction of an ecological economic system and the strengthening of the ability of environmental protection in the "Panzhihua Ecological City Construction Plan (2006~2020)". The Panzhihua Environmental Protection Bureau also put forward in the "11th Five-Year Plan" of Panzhihua Environmental Protection to vigorously develop the circular economy and create a conservation-oriented harmonious society. After 2006, ecosystem services in Panzhihua were significantly improved in the direction of synergy (Figure 8), which indicated that the coordinated ecological protection policy adopted by Panzhihua was effective.

Exploration of the Impact of Mineral Species on Ecosystem Services
For mining areas, different intensities of mining have other impacts on ecosystem services. Since the intensity of mineral extraction is difficult to judge, we explored the impact of various minerals on ecosystem services in the Panzhihua area in this study. Starting from the mining area, the ecosystem services were trade-off relationships, and ecosystem stability was undermined, which aligns with the general perception of reality. From the perspective of mineral species, the trade-off relationship was more significant in clay mines than in other mines, with the largest negative impact on ecosystem services. The study inferred that this is because clay mining requires large amounts of soil excavation, which severely destroys soil conservation, habitat quality, and carbon storage services. In contrast, water yield services result in a stronger trade-off relationship due to the pooling of water resources caused by the low-lying position of the mines. The predominant iron ore, however, has a similar trade-off synergistic relationship with the overall region and does not reflect a negative impact on ecosystem services, contrary to Xiong Jian [101], Yerner [102], etc. It is clear from the analysis that this is mainly due to the Panzhihua government's commitment to creating an ecologically protected mining area in recent years, with vanadium-titanium magnetite mines as the key construction area for both ecological and economical construction. Combined with the predicted results, there is better improvement in 2030 in terms of ecosystem service trade-offs and synergies, which indicates that the ecological construction policy of mining areas adopted in Panzhihua at the present stage has an excellent guiding effect. The Panzhihua municipal government proposed to improve the ability of coordinated and sustainable development of land and promote land development, consolidation, and reclamation in the Overall Plan of Land Use of Panzhihua (2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020). According to the prediction, the ecosystem services of Panzhihua will be better improved in 2030, which indicates that the mining ecological construction policy adopted by Panzhihua at this stage has a good guiding role. Therefore, the mining area needs to be managed in the future based on the successful policies and experiences at the current stage, with a focus on clay mining species, to grasp the development direction of the mining area and to deploy the development of the negative impact mining area to achieve the coordinated development of the ecological environment and green mining.

Uncertainty Analysis
The results in this study fill the gap in the assessment of ecosystem services in the region and, to a certain extent, provide a scientific basis for rational land use allocation and ecological protection and restoration in the future. However, at the same time, the results of this study also have some limitations. At first, the research method used in this study is limited by the accuracy of critical data such as land use classification data, which has a certain degree of uncertainty. The reliability of the study results can be further improved by using high-resolution data combined with technical means such as machine learning. Secondly, in the trade-off and synergy analysis, Spearman correlation analysis was used to explore the trade-off and synergy relationship of ecosystem services and to explore the influence brought by land use classification, but other influencing factors were not considered. In future research, we will focus on the impact of social drivers and other factors on trade-offs and synergies. Finally, due to data accuracy and timeliness in the mining impact analysis, only the effect of different mining species was investigated differently. The impact of factors such as mining intensity were not considered. In the future, we will further explore the effects of other mining factors on ecosystem services, strive to improve the analysis of mining impacts, and analyze a comprehensive mining impact analysis model to support the exploration of ecosystem services in other resource-based cities.
Panzhihua, as a typical resource-based city, is generally influenced by human activities in ecological aspects, especially in the mineral development area. Therefore, giving full play to the subjective initiative of human beings and reasonably coordinating the relationship between humans and land are essential ways to carry out the ecological restoration. At the same time, human influence is also reflected in land use, so environmental protection and restoration can be carried out by adjusting land use types. Given the situation of the Panzhihua area, this study proposes the following recommendations for coordinating the human-land relationship and ecological protection and restoration: 1.
The unsuitable arable land in the central region is supposed to be returned to forest and grass, and the mining area should be reasonably balanced between mining and replanting.

2.
Relevant authorities need to increase the construction of wetland systems, improve functional integrity, achieve synergy of multiple ecosystem services from multiple perspectives, and regulate trade-offs to build a healthier ecosystem.

3.
Most forest vegetation types in the northern part of the study area are homogeneous, and the trade-offs are more prominent. Therefore, enriching the vegetation types is necessary to improve them, enhance ecosystem stability, and promote sustainable development.

Applicability and Extension of the Model
In the study, we coupled the CA-Markov model and the InVEST model to predict future land use, evaluate the value of ecosystem services, and explore the trade-offs and synergistic relationships between ecosystem services and mining areas. This method can be extended to other spatial scales and other time scales. In view of the fact that Panzhihua is a typical mining city, we studied the impact of mineral exploitation on ecosystem services. In addition, we also adjusted the assessment unit of mining ecosystem services and proposed the assessment method based on mineral species. This provided a new way of thinking about ecosystem services and environmental assessment and governance in mining areas, depending on the type of mining.
For other ecosystems, such as grassland ecosystems, it is a good way to study grassland ecosystem services according to the functional characteristics of plants. This has been applied in the study of the relationship between plant height and edibility and ecosystem services in the grassland ecosystem of Tibet [103]. In addition, our model can also be extended to study the carbon storage service value of different forest types, or to explore the trade-offs and synergistic effects of forest ecosystem services at different north and south slope and vertical zone scales [104,105]. It should be emphasized that constructing an ecosystem services evaluation system and modifying the evaluation model combined with the actual situation of the region to improve its applicability and accuracy are the directions for further exploration in the future.

Conclusions
In this study, the CA-Markov model was used to predict the future land use of Panzhihua City, and the quantitative assessment of ecosystem service functions based on the InVEST model was conducted to explore the trade-offs and synergistic relationships among ecosystem services, and further research on the synergy and trade-offs of ecosystem services was carried out for the mining area. The results are as follows: (1) From 2002 to 2018, land use in Panzhihua has undergone dramatic changes. Cropland, water, and impervious surfaces continue to expand, and new agricultural land mainly comes from the conversion of forests and grasslands. As a result, the area of ecological lands such as forests and grasslands has decreased, and the rate of reduction has slowed down yearly. While agricultural and rural construction are vigorously developing, the conservation of ecological green areas should also be promoted. (2) The results of the ecosystem service assessment using the InVEST model showed a high degree of confidence. Water production, habitat quality, and carbon storage services have all declined, and water production services are predicted to improve in 2030. Nevertheless, the ecological condition in the economically active central part of Panzhihua is the worst, and attention needs to be paid to the adverse environmental effects of rapid urban development. (3) Synergistic relationships among ecosystem services dominate, with the most significant synergistic relationship being between habitat quality and carbon storage. Weak trade-off relationships appear between water production and soil conservation and habitat quality services, and the trade-off relationship between water production and habitat quality services is weakening. Ecosystem services in Panzhihua are gradually moving towards coordination, and ecological construction has a positive and important impact on maintaining ecosystem stability. (4) There are many trade-offs between ecosystem services in mining areas, with strong trade-offs occurring between water production and carbon storage services. Coal and iron ore mines have the most negligible impact on ecosystem services, while clay mines have the greatest impact. The effects of mining areas on ecosystem services should not be underestimated, and there is heterogeneity in the impact of different mining areas on ecosystem services. Combining ecosystem services to optimize mineral development is of great significance in achieving a win-win situation for ecology and mining. (5) A hierarchy of ecosystem service synergies and trade-offs was established. Spatially, the pattern of relationships between water yield services and soil conservation and habitat quality services showed trade-offs in the east and synergies in the west. The trade-offs between carbon storage services, water production, and habitat quality services were more significant in urban areas.
Overall, ecosystem services in Panzhihua City showed continuous improvement. The policies of mining area reclamation and promoting ecological system construction play a vital role in environmental protection. At the same time, the construction of ecological protection areas of vanadium titanium magnetite also indicates a road of harmonious development for mineral mining and ecological construction. It is particularly important to actively promote the policy construction of ecological restoration and protection.
Panzhihua is a resource-based city with rich and diverse mineral resources. We adjusted the model and proposed an evaluation method based on multiple minerals. Our study on ecosystem service trade-offs and synergies for nine key minerals, including coal and iron, could provide a reference and demonstration for other mining cities.