Spatiotemporal Variation in Ecosystem Services and Their Drivers among Different Landscape Heterogeneity Units and Terrain Gradients in the Southern Hill and Mountain Belt, China

: Understanding the spatiotemporal heterogeneity of ecosystem services (ESs) and their drivers in mountainous areas is important for sustainable ecosystem management. However, the effective construction of landscape heterogeneous units (LHUs) to reﬂect the spatial characteristics of ESs remains to be studied. The southern hill and mountain belt (SHMB) is a typical mountainous region in China, with undulating terrain and obvious spatial heterogeneity of ESs, and was selected as the study area. In this study, we used the fuzzy k-means (FKM) algorithm to establish LHUs. Three major ESs (water yield, net primary productivity (NPP), and soil conservation) in 2000 and 2015 were quantiﬁed using the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model and Carnegie Ames-Stanford approach (CASA) model. Then, we explored the spatial variation in ESs along terrain gradients and LHUs. Correlation analysis was used to analyze the driving factors of ESs in each terrain region and LHU. The results showed that altitude and terrain niche increased along LHUs. Water yield and soil conservation increased from 696.86 mm and 3920.19 t/km 2 to 1061.12 mm and 5117.90 t/km 2 , respectively, while NPP decreased from 666.95 gC/m 2 to 648.86 gC/m 2 . The ESs in different LHUs differed greatly. ESs increased ﬁrst and then decreased along LHUs in 2000. In 2015, water yield decreased along LHUs, while NPP and soil conservation showed a ﬂuctuating trend. Water yield was mainly affected by precipitation, temperature and NDVI were the main drivers of NPP, and soil conservation was greatly affected by precipitation and slope. The driving factors of the same ES were different in different terrain areas and LHUs. The variation and driving factors of ESs in LHUs were similar to some terrain gradients. To some extent, LHUs can represent multiple terrain features. This study can provide important support for mountain ecosystem zoning management and decision-making.


Introduction
Ecosystem services (ESs), which are coupled with natural and social systems, directly and indirectly influence human survival and development [1,2]. ESs link ecosystems with economic activities and have become a hot topic of research in recent years. Additionally, ESs are closely linked to human well-being and can form important components of sustainable management policies [3][4][5]. Decision-makers seek to maximize ESs through effective management [6]. Therefore, identifying the spatiotemporal variation characteristics of ESs is critical for designing interventions. Many scholars have explored the spatiotemporal differences in ESs in different regions and different environmental gradients [7][8][9]. Landscape heterogeneous units (LHUs) are areas with similar spatial heterogeneity and form the basis of landscape management [10][11][12]. However, few studies have considered ESs across different LHUs. Although some researchers have discussed the spatial heterogeneity of ESs among different terrain gradients [7,13,14], this research is not sufficient.
As one of the most important providers of nature's contribution to humans, mountain ecosystems are fragile and subject to multiple drivers of change [15]. Biophysical conditions, including temperature, precipitation, evapotranspiration, soil type, and ultimately vegetation type, can vary dramatically in mountain systems [16]. Therefore, the spatial heterogeneity of ESs is high [17]. However, we still have a limited understanding of how to construct landscape units to reflect spatial heterogeneity and how provisioning patterns of multiple ESs emerge and change in space and time along the landscape and terrain gradients in mountain systems.
Identifying the mechanisms underlying the drivers of ESs can provide important information for the development of targeted ecological management policies [18][19][20][21][22]. Therefore, we can improve multiple ESs by manipulating the drivers of ESs. Most previous studies have focused on the spatial distribution, mapping, and relationships of ESs at a single point in time [23,24]. Recently, an increasing number of researchers have investigated the complex relationships among ESs and how their relationships might change over longer timespans [25][26][27][28][29] by analyzing land cover/land use (LULC) changes [30][31][32], performing scenario analysis [33,34] or assessing ES bundle changes [35,36]. The overall research on ESs is in-depth and specific, with a variety of perspectives. The relationships, changes and spatial heterogeneity of ESs often depend on driving factors [13,14]. Some scholars have explored the spatial heterogeneity of driving factors [37][38][39]. However, how different drivers change across terrain and the heterogeneity of landscape units has been less considered, and few studies coupling socioecological factors in the context of heterogeneous landscapes have explored the driving mechanisms behind ESs [39]. Therefore, understanding these drivers in different LHUs is urgently needed.
Ecosystem responses to direct drivers vary along spatial and temporal scales and affect ecosystem goods and services at different scales [18,40,41]. The diverse landform types, substantial topographic relief, and strong spatial heterogeneity make the mountain ESs spatially diverse [7,8,13,14]. For instance, Sun explored the trade-offs and synergies among multiple hydrological ecosystem services (HESs) in different topographic basins, and the results demonstrated that the intercorrelations among HESs varied significantly between the western plain and eastern mountain areas [13]. Wang determined that the ESs showed high spatial heterogeneity in different geomorphological types due to the differences in the intrinsic characteristics of each geomorphological type in the Beijing-Tianjin-Hebei region of China [14]. Additionally, Gomes showed the spatiotemporal variation in the contrasting responses of ESs to different biophysical conditions in the three altitude zones along an altitudinal gradient in Brazil [17]. Therefore, a major research challenge is to understand the mechanisms that form and change ESs at different landscape units. To this end, it is necessary to construct multiple LHUs and analyze the spatial and temporal changes in ESs to provide effective information for sustainable ecosystem management.
The fuzzy k-means (FKM) algorithm is a widely used clustering technique [42], and it is regarded as an effective method that can produce results with high within-cluster homogeneity and between-cluster heterogeneity [43]. In the field of ES research, FKM has been used to cluster ES bundles [35,[44][45][46], rural-urban classification [47,48], and major vegetation clustering [8]. Few studies have considered obtaining multi-LHUs using the FKM, which maximizes the cluster heterogeneity of the landscape. Therefore, in this study, we selected the southern hill and mountain belt (SHMB), which is an important part of the "two shelters and three belts" in China, as the study area. First, based on 5 terrain parameters, we constructed the LHUs using the KFM algorithm. Second, we estimated 3 key changes in ESs across the different LHUs and terrain regions in 2000 and 2015. Third, we used correlation analysis to analyze influencing factors of ESs in each LHU and terrain region. Finally, we proposed ecological management suggestions for different LHUs. The results can provide theoretical guidance and decision support for landscape management in heterogeneous mountainous areas.

Study Area
The southern hill and mountain belt comprises 114 counties with an area of 288,655 km 2 , which includes 6 provinces: Yunnan Province, Guizhou Province, Hunan Province, Guangxi Province, Jiangxi Province, and Guangdong Province (Figure 1). In this area, the elevation increases from the southeast to northwest, with the highest elevation of 2861 m. The study area has a subtropical monsoon climate with average annual rainfall ranging from 1400 mm to 1800 mm. The region has a variety of vegetation types, and forests account for more than 70% of the land cover.
(3) The digital elevation model (DEM) was downloaded from the Resources and Environmental Science, Chinese Academy of Sciences (http://www.resdc.cn/ accessed on 24 September 2019), with a spatial scale of 1000 m × 1000 m. The terrain parameters such as altitude, slope, aspect, landform relief, and terrain niche were derived from the DEM using SAGA software (http://www.saga-gis.org/en/index.html accessed on 24 September 2019). These parameters were used for cluster analysis to obtain the LHUs.

Water Yield
The water yield of different landscape types was calculated based on the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) water yield module based on the principle of water balance [49]: where Y xj is the water yield of pixel x for LULC type j, AET xj is the actual evapotranspiration of pixel x for LULC type j, and P x is the annual precipitation of pixel x. The input parameters and calculation methods of the water yield model refer to the studies of Fang and Bai [8,30].

Soil Conservation
Soil conservation refers to the ability of an ecosystem to control soil erosion and retain sediment. Soil conservation was calculated based on the universal soil loss equation using the sediment delivery ratio (SDR) module of the InVEST model [30,50]. The formula is as follows: where A is the soil conservation capacity, R is the precipitation erosivity, K is the soil erodibility, L is the slope length factor, S is the slope gradient factor, C is the vegetation management factor, and P is the practice factor. L and S were calculated by the InVEST model. The biophysical table and related parameters, including R, K, C, and P, were derived from different studies [51][52][53]. Other parameters, such as L and S, were calculated by the InVEST model.

Net Primary Productivity
In this study, we used the Carnegie Ames-Stanford approach (CASA) model based on the light energy utilization principle to estimate NPP [54,55] as follows: where APAR(x, t) is the effective photosynthetic radiation absorbed by pixel x in month t, and ε(x, t) is the actual utilization rate of light energy.
APAR(x, t) = SOL(x, t) × FRAR(x, t) × 0.5 (4) where SOL(x, t) is the total solar radiation of pixel x of month t and FRAR(x, t) is the absorption ratio of vegetation to incident photosynthetically active radiation.
where T ε1 (x, t) is the inhibitory effect of low temperature on light utilization, T ε2 (x, t) is the inhibitory effect of high temperature on light utilization, and W ε (x, t) is the stress coefficient affected by water. The parameter input of the CASA model used data from Zhu [55].

Terrain Analysis
On the basis of relevant studies, we selected altitude, slope, aspect, landform relief, and terrain niche to analyze the characteristics of ecosystem services. In this study, we divided the study area into five regions (each approximately 20% of the total area) according to terrain parameters (Table 1). Random samples of terrain gradients were then selected for analysis of variance (ANOVA) to reveal differences in ESs.

Landscape Heterogeneity Units
The fuzzy k-means (FKM) algorithm calculates the membership of objects belonging to specific clusters, and the degree of similarity among the objects can be easily obtained. Therefore, the FKM algorithm can be used to establish spatial LHUs and zoom in on landscape differences.
The function of the FKM algorithm is defined by Equations (6)- (8): where J m is the objective function, n is the number of objects in the database, c denotes the number of clusters, m is the fuzzy exponent, u ij is the degree of membership of object X i in cluster j, and k j is the d-dimensional center of the cluster (Dipartimento di Elettronicae Informazione). By updating the values of the membership u ij and the cluster center k j , fuzzy partitioning iteratively optimizes the objective function. The routine stops when the objective function J m converges to a minimum. The fuzzy exponent (m), which defines the fuzziness among each cluster, and the number of clusters (k), which indicates the total number of clusters, are major parameters used in FKM. In this study, we chose five clusters (k = 5), which were equal to the five terrain gradients, and m was equal to 1.2 for cluster analysis.

LULC in Different Terrain Gradients and LHUs
From 2000 to 2015, LULC in SHMB was consistent, with a change of less than 1% [53]. Therefore, we only analyzed the spatial variation of LULC. Table 2 shows the composition of the LULC for each terrain gradient and LHU. FL occupied the largest proportion (41.41%) of the SHMB. FL comprised the largest proportion of the LULC in all terrain zones. SL, OC, GL, PF, and DF were also the main LULC types in the SHMB (13.91%, 13.31%, 12.33%, 8.97%, and 7.92%, respectively) ( Figure 2). The proportion of SL in the high-altitude zones was higher than that in the low-altitude zones, while the proportion of SL was the opposite with respect to aspect gradients. SL occupied a high proportion at slopes of 1.16 • -2.33 • (3.46%), landform reliefs of 78-139 m (3.60%), and terrain niches of 1.10-1.41 (3.16%). OC decreased with the increase in aspect. OC occupied a large area at altitudes of 680-1028 m (3.08%), slopes of 2.33 • -3.84 • (2.98%), landform reliefs of 139-208 m (2.89%), and terrain niches of 1.41-1.74 (3.18%). GL tended to increase almost with the terrain gradients except in aspect zones. GL comprised a large proportion of 134.89 • -212.75 • (2.89%). The distributions of PF and DF were similar, and more areas of these LULC types were distributed in the low-terrain zones than in the high-terrain zones.  The spatial distributions of terrain gradients and LHUs are shown in Figure 3. Unlike terrain zones, the area of LHUs was different. The areas of the first to fifth LHUs were 95,319, 78,587, 59,266, 31,537, and 24,081 km 2 , respectively. From the first to the fifth LHUs, the altitude and terrain niche showed increasing trends. However, the aspect of each LHU exhibited no obvious change. Slope and landform relief showed a trend of first increasing and then decreasing from the first to the fifth LHU and reached their peak at the fourth LHU. FL dominated the main components of the first through third LHU (12.93%, 15.03%, and 9.07%, respectively) and occupied a large area in the fourth LHU. SL and GL first increased and then decreased along the LHUs, reaching their peaks in the second (3.47%) and third (3.10%) LHUs, respectively. OC was more distributed in the first three LHUs (3.16-3.97%). PF and DF tended to be concentrated in the first LHU (5.81% and 3.23%, respectively).

Spatiotemporal Variation in ESs from 2000 to 2015
ESs showed different spatial patterns over time except for soil conservation (Figure 4), which was in line with the results of Yin (2019) [56]. The average water yield in the SHMB in 2000 was 696.86 mm, showing a spatial distribution of high water yield in the east and low water yield in the west. In 2015, the average water yield of the SHMB increased to 1061.12 mm, and the high value of water yield was concentrated in the central part of the SHMB, while the low values of water yield were located in the western and eastern parts of the SHMB. The water yield of the SHMB from 2000 to 2015 showed an increasing trend overall, but there were obvious regional differences. From 2000 to 2015, the decrease in water yield mainly occurred in small parts of the southwestern and southeastern SHMB, while the increasing trend of water production occurred in most areas, especially in the central region of the SHMB. In particular, 7.36% of the area exhibited decreased water yield, while 91.37% of the area showed increased water yield, in the SHMB.
The  In contrast to water yield and NPP, the spatial distribution of soil conservation in the SHMB did not change significantly from 2000 to 2015. The average soil conservation in the SHMB in 2000 was 3920.19 t/km 2 , which increased to 5117.90 t/km 2 in 2015. Spatially, the soil conservation in the eastern part of the SHMB was slightly higher than that in the western part of the SHMB. Specifically, soil conservation increased in 81.38% of the SHMB and decreased in 8.59% of the SHMB. In addition, the decrease in soil conservation was mainly concentrated in the eastern part of the SHMB. Overall, ESs showed different trends and had obvious regional spatial differences.

Changes in ESs in the Terrain Zones and LHUs
The mean ES values in each terrain zone and LHU are shown in Figure 5. Differences in the mean levels of ESs were assessed using ANOVA to reveal differences in ESs (p < 0.05) across terrain zones and LHUs. The water yield in 2000 in the three lowest altitude zones was similar (739. 12  In particular, for the same terrain zone or LHU, the water yield in 2000 was always lower than that in 2015. In general, the water yield was not as good in high-terrain zones as in low-terrain zones, while the soil conservation was the opposite. High NPP was mainly distributed in medium-terrain zones. For LHUs, the third and fourth LHUs had high water yields, the first and second LHUs had high soil conservation, and the second and fourth LHUs had high NPP.

Spatial Patterns and Drivers of ESs in Terrain Zones and LHUs
Understanding the spatiotemporal variations and drivers of ESs can provide important information for ecological management [17,30,50]. Many scholars have studied ESs in mountainous areas, but few studies have evaluated ESs in complex terrains to explore the differences in ESs among different terrain regions [7,13,57]. Our results showed that there were great differences in water yield with altitude, terrain niche, and LHUs but small differences with slope, aspect, and landform relief. In contrast, NPP and soil conservation had obvious differences in terrain regions and LHUs except for aspect. In addition, different LHUs had terrain gradient effects except for aspect. Therefore, terrain is an important basis for the distribution of ESs, but aspect cannot be considered when analyzing the terrain gradient effect of ESs in the SHMB. Effective information on the terrain gradient effect of ESs can provide a scientific reference for the coordinated development of social ecology and the optimization of ESs.
The spatial patterns of ESs in mountainous areas are significantly different, and ecological management should be adapted to local conditions [8,58]. Significant differences in climate, terrain, LULC, and other socioecological factors in the SHMB resulted in different spatial patterns of ESs [53]. In particular, we selected precipitation, temperature, Normalized Vegetation Index (NDVI), and slope as influencing factors and used correlation analysis to reveal the relationship between ESs according to relevant literature and influencing factors in different terrains and LHUs (Table S1). The distributions of precipitation and water yield were relatively similar. In addition, correlation analysis showed that precipitation was always the main factor impacting water yield in all terrain regions and LHUs (r > 0.361, p < 0.01). Moreover, some studies also pointed out that there is a strong positive correlation between the changes in water yield and precipitation, and the spatial distribution of water yield is highly consistent with precipitation [59][60][61]. Thus, the increased precipitation led to the significant increase in water yield. Compared with 2000, the spatial distribution of precipitation in the SHMB in 2015 changed significantly, leading to a great change in the spatial distribution of water yield in 2015. The composition of LULC also affected ESs [62,63]. The FL in low-terrain areas was less than that in high-terrain areas, and the water yield of FL was lower, which also led to the terrain gradient change in water yield. There was a significant negative correlation between water yield (r < −0.180, p < 0.01) and temperature (r < −0.288, p < 0.01) at different elevations and LHUs in 2015, which indicated that temperature had a negative effect on water yield. Water yield was positively correlated with the NDVI (r > 0.090, p < 0.01) in all terrain regions in 2000. This result may have occurred because vegetation with higher coverage has more canopy interception and ground litter water capacity [64]. The spatial distribution of NPP was positively correlated with temperature (r > 0.106, p < 0.01), indicating the contribution of temperature to NPP and NPP was more sensitive to temperature. NPP was positively correlated with the NDVI (r > 0.102, p < 0.01) in many regions, which was obvious because the higher the vegetation coverage was, the stronger the carbon sequestration ability was. Temperature and NDVI were the main factors affecting the spatiotemporal distribution of NPP, which was consistent with other studies [65,66]. Changes in temperature and spatial distribution of NDVI led to changes in NPP. Terrain factors also influence the distribution of vegetation zones by affecting the vertical zoning of climate, which also affects NPP. The precipitation in SHMB is abundant, and the area with relatively steep terrain is susceptible to soil erosion and the resulting destruction of the regional vegetation. In addition, Fu pointed out that the central part of the SHMB is suffering from ecosystem degradation caused by deforestation, which also confirms the decline in vegetation coverage, resulting in a decrease in NPP [51]. This effect explains the decline in NPP in the SHMB. The relationship between precipitation and NPP was unstable, with both positive and negative correlations, indicating that rainfall had an inconsistent influence on NPP. The significant positive correlation between soil conservation and slope highlighted the effect of slope on soil conservation in mountainous areas [28]. In addition, the NDVI and precipitation were high in areas with high soil conservation, which was consistent with previous studies [64]. Precipitation is also a major factor affecting soil conservation. Changes in precipitation distribution led to changes in soil conservation distribution. An increase in precipitation directly leads to an increase in rainfall erosivity and a corresponding increase in soil conservation. Precipitation, temperature, NDVI, and slope were spatially heterogeneous, leading to spatial differences in ESs. Comparing the driving factors of ESs on terrain gradients and LHUs revealed that they were basically the same. Therefore, it is feasible to use LHUs instead of terrain gradients to analyze the distribution and driving factors of ESs. Understanding the drivers in different LHUs is critical for policy development and management.
Interestingly, the relationship between the same factor and water yield was different in different areas of the same terrain type, as shown by the relationship between NPP and soil conservation. For example, the relationship between temperature and water yield presented a positive correlation (r = 0.130, p < 0.01), no correlation (r = 0.051, p > 0.01), and a negative correlation (r < −0.247, p < 0.01) from low altitude to high altitude in 2000. Similar findings have been found in other topographic regions and in LHUs, suggesting that driving factors change with terrain gradient, LHUs, and time, underscoring the importance of considering driving factors in different terrain areas and LHUs [13,14].
ESs and their drivers have terrain gradient effects; however, many current ecological conservation, restoration and management policies tend to take a macro view and ignore the spatial heterogeneity caused by terrain factors. However, there are many terrain parameters, and some terrain indices are not applicable to a region. In addition, considering multiple terrain parameters at the same time presents considerable challenges for decision-makers. The concept of LHUs proposed in this study can simplify the terrain parameters considered in policy, and the results showed that ESs also had significant spatial variation in LHUs and even greater variation in water yield. Moreover, decision-makers can select specific terrain indices to obtain LHUs according to actual needs. Therefore, it is of great value to carry out effective ecosystem management to understand the impacts of terrain factors on the spatial differences and driving factors of ESs in relation to LHUs.

Ecological Management Measures in Different Regions
As an ecological barrier area in China, the SHMB is also the headwaters of the Yangtze River and the Pearl River, which are of great importance for ensuring ecological security and promoting sustainable development in southern China [51]. The SHMB has implemented the national major ecosystem protection and restoration project master plan, natural forest protection program, and Grain for Green Program to improve the natural ecosystems [51,53]. However, insufficient attention has been given to the changing terrain, making these programs difficult to implement; thus, the SHMB is now facing ecosystem degradation from deforestation [51]. Thus, combined with the terrain effects of LULC, ecosystem services and their drivers, the following targeted ecological management measures in different LHUs are proposed.
The low elevation of the first LHU and the abundant water supply provide favorable conditions for agricultural production, thus driving a large and concentrated population. The first LHU is the main agricultural supply region in the SHMB, but FL occupies the largest proportion of the area, and FL has high ESs. Therefore, we should prioritize the protection of existing vegetation, prohibit commercial logging and deforestation for land reclamation, and develop organic agriculture to increase crop yields [8,67]. This region is an important area for economic development and is also the most densely populated area in the SHMB. Thus, ecologically sustainable construction is necessary for the common needs of ecological protection and economic development to seek sustainable ecosystem management. Commercial plants can be incorporated into mixed agriculture to reduce ecological impact, and forest parks can be established to promote forest tourism [8,68]. In addition, against the background of global warming, the trend of extreme precipitation events has generally emerged, and extreme precipitation events will have many negative impacts on the ecological environment [69]. As the region with the most precipitation in China, it is necessary to strengthen the monitoring of precipitation in the SHMB in the future to provide a scientific basis for regional ecosystem management to cope with future climate change.
In the second and third LHUs, the quality and stability of forest ecosystems should be enhanced. The area is dominated by natural vegetation. We should comprehensively protect evergreen broad-leaved forests and other natural vegetation, improve the quality of forests in a scientific and precise manner, nurture young and middle-aged forests, and restore degraded forests. We should also establish ecological reserves, identify and connect ecological corridors, and improve biodiversity protection networks [70]. Additionally, we can steadily carry out the natural forest protection program, the Grain for Green Program, and the task of controlling soil erosion. The government can appropriately encourage the development of understory planting and mixed forest to improve the use of regional land resources [71]. Agroforestry systems can consider both natural and socioeconomic objectives, improve the local microclimate, control soil erosion, improve ESs and increase farmers' incomes [72,73]. Demonstration zones for agroforestry, forest tending, and restoration should be established and promoted. Ecological supervision should be enhanced, and the boundaries of natural vegetation should be effectively monitored to prevent their destruction.
The fourth and fifth LHUs are the most important ecological restoration areas in the SHMB. This area is a typical karst area and a relatively poor area. The comprehensive control of soil erosion and rocky desertification should be promoted, and the ecological restoration of mines and comprehensive land improvement should be carried out. The restoration of FL and GL and the Grain for Green Program need to be intensified [74].
Farmers should be encouraged to develop animal husbandry, plant high-quality forage grass on sloped farmland, and build a complex agro-pastoral complex ecosystem. Complex agro-pastoral ecosystems can balance ecological restoration with the development of characteristic industries and realize harmony between humans and nature in fragile geological environments [75]. It is necessary to accelerate the construction of spatial supervision and performance evaluation systems to promote the effective management of different types of ecosystem. A demonstration area of poverty alleviation through science and technology should be established to promote overall poverty alleviation and rural revitalization.

Limitations and Future Directions
We compared the model outputs of water yield, NPP, and soil conservation with those of other studies in the SHMB. Yin (2019) estimated the perennial average water yield, NPP, and soil conservation in the SHMB to be 721.4 mm, 672.4 gC/m 2 , and 4287.4 t/km 2 , respectively, which were similar to the results of this study in 2000 [56]. The precipitation in 2015 was significantly higher than that in 2000, and the water yield and soil conservation were greatly affected by precipitation; thus, the water yield and soil conservation in 2015 were relatively large. Due to the different time ranges and regional scales of different studies, there are some differences in the results, but the calculation results in this study were broadly within the range reported in the study of Yin (2019). Since LULC data are used in ES assessments, the accuracy of LULC can lead to bias in ESs. In the future, researchers should improve the classification accuracy of LULC to reduce errors. The InVEST water yield module ignores the interaction between surface water and groundwater and the upstream water recharge, which may lead to errors in water yield calculations. In addition, the lack of long-term field observation data in this study would affect the accuracy of the results. In future research, the monitoring and collection of field data should be strengthened, and the model parameters should be adjusted to ensure the reliability of the evaluation results. In the SDR module of the InVEST model, the same LULC type has the same C value and p value, while the same LULC type is different in space, which would lead to errors in estimations of soil conservation. In the CASA model, the maximum light energy utilization rate used in this study is a universal value, while the maximum light energy rate of different LULC types is significantly different, which will bring errors to the estimation of NPP. Therefore, in future studies, it is necessary to conduct in-depth discussions on the maximum optical energy rates of different LULCs to determine the appropriate maximum optical energy rate parameters. LHUs are defined by the FKM algorithm based on terrain parameters, so the deviation of LHUs mainly comes from the accuracy of DEM and the selection of terrain parameters. In the future, stakeholder analysis can be used to select the appropriate terrain parameters and use high-precision DEM to construct LHUs.
In this study, only three key ESs in the SHMB were considered, and multiple ESs could be reasonably selected according to the characteristics of the study area in the future. This study considered only the ESs in 2000 and 2015 at two time nodes; there was a lack of time and space dynamic analysis of time series for a long continuous period, and the ESs fluctuated with time. Future studies should involve a long time series data analysis of ES changes in LHUs to reveal the general trends, as longer periods are needed for an ecological system to exhibit reductions in abnormal factors-and due to the time lag effect caused by errors; this longer-term perspective would improve the reliability of the analysis [76]. The results provide insights into the terrain gradient effect of ESs and their influencing factors, but the driving mechanisms remain to be explored to address future challenges. The driving factors considered in this study are relatively one-sided, and in the future, climate factors, terrain factors, socioeconomic factors, and land-use structure factors can be comprehensively considered, and their temporal and spatial differentiation can be estimated. ESs often have complex trade-offs/synergies, and whether these relationships change with terrain gradients remains to be studied [13]. Future studies may explore whether the tradeoff/synergy of ESs changes with LHUs. Ecological management policies that do not consider the interactions of ESs may result in an increase in one ES but a decrease in another. It is necessary for future research to predict ES interactions and terrain gradient changes under future land-use management scenarios and determine the key factors affecting ESs to ensure regional ecological environment security and achieve a win-win situation for ecological and social development.

Conclusions
A comprehensive consideration of the impact of terrain factors on the spatial differences in ESs and their drivers can provide useful information for policy formulation. The spatial heterogeneity of ESs may be further understood through the construction of LHUs. In this study, LHUs were constructed based on terrain parameters according to the FKM algorithm, and ESs were quantified through InVEST and CASA models. Then, spatiotemporal variations of ESs along terrain gradients and LHUs were analyzed. Correlation analysis was used to analyze the driving factors of each terrain and LHU. The results showed that the spatial distribution of ESs was heterogeneous, and the spatial distributions of water yield and NPP changed greatly over time. The variation of ESs with LHUs is similar to that of some terrain parameters, indicating that LHUs can reflect terrain characteristics. Compared with other terrain parameters, the effect of aspect on the difference in ESs was small. Suitable terrain parameters should be selected according to the characteristics of different regions. Correlation analysis showed that precipitation and temperature were the main factors affecting water yield, temperature and the NDVI affected NPP, and precipitation and slope were the main drivers of soil conservation. Changes in driving factors have obvious temporal and spatial differences. Identification of the driving factors of different LHUs can provide a reference for ecosystem protection and restoration in different regions. The similarity of spatiotemporal differentiation and drivers of ESs between LHUs and terrain gradients indicated that LHUs could be used to represent multiple terrain parameters. Thus, LHUs can help decision-makers to consider the spatial heterogeneity of ESs when making decisions. According to the characteristics of each LHU, we proposed different ecological management measures. The method of using FKM to establish LHUs in this study can be applied to other regions and incorporated into the formulation of ecological zoning management policies in mountainous areas. Therefore, the relationship between landscape heterogeneous units and ecosystem management should be considered in future studies.

Supplementary Materials:
The following is available online at https://www.mdpi.com/article/10 .3390/rs13071375/s1, Table S1: Correlation between the ESs and influencing factors along the terrain gradients and LHUs.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.