The Hengduan Mountain region (24°39′N-33°34′N, 96°58′E-104°27′E) is located in Southwest China and has a total area of 449,748 km² that mainly includes the eastern part of the Tibet Autonomous Region, western Sichuan Province, and northwestern Yunnan Province. The elevation ranges from 306 to 7143 m and tends to increase from southeast to northwest. This region also contains various geomorphological zones because of the complex terrain (Figure 1). The study area is an upstream area of the major rivers in China and Southeast Asia, with developed river systems and numerous tributaries. The climatic zoning of this region includes mainly the Qamdo and Bomê-western Sichuan regions belonging to the plateau climate zone, and the Sichuan, northern Yunnan, Jinsha river-Chuxiong, Yuxi regions that belong to the mid-subtropical climate zone (Figure 1). The vegetation also presents an obvious vertical heterogeneity due to the distinct differences in climatic conditions, which vary from the arid valley shrubs at low elevations to the sparse alpine vegetation found at high elevations (Zhang et al., 1997). Additionally, the vertical zonation of the vegetation varies based on the horizontal zone. The soil types cover four major soil series: red, brown, cinnamon, and alpine soil. The spatial differentiation of vegetation and soil and the various combinations of the two result in changes in the water yield service across the study space by affecting evapotranspiration and other water cycle processes, finally leading to the spatial differentiation of water yield.

According to the water yield model in the InVEST software and geographical detector method, the data required in this study mainly included climate data, a digital elevation model (DEM), soil data, land use data, a normalized difference vegetation index (NDVI), hydrological data, geomorphological, and climatic types data. Climatic data such as daily temperature, precipitation, and wind speed of the 42 meteorological stations—located in and around the Hengduan Mountain region, were obtained from the National Meteorological Information Center of China (http://data.cma.cn/); Raster-formatted meteorological data were obtained via Anusplin interpolation method; DEM (mainly used for the extraction of watersheds and sub-watersheds), geomorphological and climatic types data from the resource and environment data cloud platform (http://www.resdc.cn/). Soil data included soil type, depth, texture, and organic matter content. The land use data used in this study was obtained viainterpretation of 30 m resolution remote sensing images. We divided the land use types into 22 categories (Figure 1) as per the classification system of Liu et al. (2010) with an accuracy of 91.04%, to meet the requirements of our study. Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data (MOD13A3: 1 km resolution monthly data) was obtained from the Land Processes Distributed Active Archive Center (https://lpdaac. usgs.gov/). Based on the land geomorphology of Zhou et al. (2009), we merged the types according to the relief degree of land surface and formed seven types of terrains namely, plains, platforms, hills, small relief mountains, medium relief mountains, large relief mountains, and extreme relief mountains (Figure 1). The runoff data during 2006-2015 for the 18 hydrological stations were obtained from the hydrologic data yearbook of the People’s Republic of China and used for the verification of the model-simulated results.

2.3.1 Framework

The present study focused on the spatial heterogeneity of water yield service in the Hengduan Mountain region and investigated the driving factors underlying this heterogeneity. Based on the above objectives, we constructed the following research framework (Figure 2): (1) Based on the water yield model in the InVEST software, we evaluated the water yield service and validated the results. (2) The spatial differences in the degree of relief and the climate are the basis for the heterogeneity of various ecological and environmental elements in mountain areas. We explored the main influencing factors for the spatial heterogeneity of water yield in different geomorphological and climatic zones by using the geographical detector method. (3) We identified the contribution of single factors to the water yield heterogeneity with a ‘factor detector’ and an ‘ecological detector’, determined the effects of the interactions of different factors with an ‘interaction detector’, and identified the areas of high value water yield service with a ‘risk detector’ .

2.3.2 Water yield service assessment and mapping

The water yield module in the InVEST software is based on the Budyko curve (Budyko, 1974):

where Y_xj and AET_xj are the average water yield and annual actual evapotranspiration on the pixel x for the land use typej, respectively, P_x is the precipitation on the pixelx. AET_xj/P_x is the ratio of the actual evapotranspiration and precipitation which is based on the development of Budyko curve by Zhang et al. (2001):

where R_xj represents the Budyko dryness for the pixel x and land use type j; ${\omega _x}$ is the ratio of plant accessible water storage and the expected precipitation; Z is Zhang coefficient (Zhang et al., 2001); k_ij is vegetation evapotranspiration coefficient; $PAW{C_x}$ is plant available water content (Zhou et al., 2005).

2.3.3 Identifying influencing factors for the spatial heterogeneity of water yield

The geographical detector method is a set of statistical methods that detect the relationship between an event and potential risk factors (Wang and Xu, 2017). If the spatial distributions of the independent variable and the dependent variable are similar, then the independent variable has an important influence on the dependent variable (Wang et al., 2010; Wang and Hu, 2012). Greater levels of similarity imply a greater influence. In this study, we set the water yield service as the dependent variable and identified the main factors for its heterogeneity by using the geographical detector method in different geomorphological and climatic zones. In terms of influencing factors, this study focused on both natural and human activities. The natural factors involve climate, topography, soil, and vegetation, and human activity factors mainly refer to those related to land use. Finally, 15 factors were selected as independent variables (Table 1). In addition to the land use type, the factors related to the composition and structure of land use had been considered as land-use factors expressed by landscape metrics. Specifically, at the landscape level, Shannon’s diversity index (SHDI), Contagion Index (CONTAG), and Effective mesh size (MESH) were selected to represent the diversity, concentration, and fragmentation, respectively. Class level metrics of the Percentage of Landscape (PLAND) for cultivated land, forests, and grassland were used to present the land use constitution. All the metrics were calculated by the ‘moving window’ function in the landscape pattern analysis model Fragstats (Mcgarigal and Marks, 1995). A total of 16,925 points were generated at 5 km intervals for the spatial correlation of multiple factors.

In general, the geographical detector contains four formulas: the factor detector, the interaction detector, the risk detector, and the ecological detector. The factor detector is used to detect the spatial differentiation of water yield service and to determine the proportion of the spatial distribution of water yield that can be explained by different factors. It is obtained by comparing the sum of the variance of the subareas in a region and the variance of the total region, which can be measured by the q value (Wang et al., 2010):

where h indicates the stratified status of water yield or climate, topography, soil, vegetation, land use, and other factors, and there were a total of L strata. The input data required in geographical detector must be categorical data, therefore, we divided the continuous data into strata data. With reference to the data discretization method and the study experience proposed by Wang and Xu (2017), the slope data was divided into 6 levels (≤5°, 5°-8°, 8°-15°, 15°-25°, 25°-35°, >35°), NDVI was divided into 7 levels according to the natural breaks, and the other factors were divided into 9 levels with the same method. All the influencing factors used unitive discretization standards across different geomorphological and climatic zones. SSW andSST represent the sum of the variance of the subareas and the variance of the total region, respectively. N is the number of units in the whole region, and N_h is the number of units inh stratum; $\sigma _h^2$ is the variances of water yield in theh-stratum, and ${\sigma ^2}$ is the variances in the whole region. The q value indicates that influencing factors could explain 100×q% of the spatial distribution of water yield. A larger value means a stronger explanatory power whose significance can be detected by the geographical detector software (Wang and Xu, 2017).The ecological detector is used to determine if there is a significant difference in the contribution of different influencing factors to the water yield spatial heterogeneity. The F statistic was usually used to detect when the q values of two factors were not significantly different:

In the risk detector, the t test was used to explore significant differences in water yield across different strata:

where ${\overline Y _h}$ indicates the mean value in the h stratum, and Var represents variance.The contribution of the interaction of two influencing factors to water yield can be calculated using the interaction detector module. If q (X1∩X2) < Min (q(X1), q(X2)), the interaction presents nonlinear weakening; if Min (q(X1), q(X2)) < q(X1∩X2) < Max(q(X1),q(X2)), the interaction is single factor nonlinear weakening; if q(X1∩X2) > Max(q(X1),q(X2)), the interaction is two factors enhancement; if q (X1∩X2)=q(X1) + q(X2), the two factors are independent; if q(X1∩X2) >q(X1) + q(X2), the relationship is nonlinear enhancement.

The water yield depth of the Hengduan Mountain region ranged from 0 to 1619 mm in the year 2010, and the average depth of the water yield of the whole region was 420 mm. Based on the DEM data, we rebuilt the corresponding sub-watersheds by taking the hydrological stations as watershed outlets through ArcHydro. Further, we compared the observation runoff data with the simulated data of each sub-watershed for verification, and the results showed that there was a high correlation (the coefficient reaches as high as 0.9947). This indicates that the simulated results can be used to analyze the spatial heterogeneity of water yield service—in which spatial differences exist. The high water yield value areas were mainly concentrated in the Three Parallel Rivers Region (TPRR) in the western part of the Hengduan Mountain region and in areas of Yuexi, Meigu, and Zhaojue counties located at the southwest edge of Sichuan Basin. The Chaya, Gongjue, and Mangkang county areas lying in the northwestern part belong to the areas of low water yield value (Figure 3).

For the Hengduan Mountain region, precipitation (58.9%) and actual evapotranspiration (26.7%) had the strongest explanatory power for water yield, which were significantly higher than those of other factors. The ecological detector showed that the effects of precipitation and actual evapotranspiration on water yield were significantly different (Figure 4). For factors with q value between 10% and 15%, the explanatory degree decreases in following order: land use type (14.7%), solar radiation (14.5%), elevation (12.6%), and soil type (11.08%). Significant differences existed in the effects of most pairwise factors on water yield except for land use type with solar radiation and elevation with soil type. Among the factors with an explanatory degree of less than 10%, NDVI and annual mean temperature had higher explanatory degree, and there were no significant differences between the two. We classified the factors and found that the q values of climatic factors were the highest and that significant differences existed among these factors. For land use related factors, the contribution of land use type to water yield was much higher than the contribution of vegetation growth (NDVI), land use composition, or structure. In terms of topography related factors, elevation seemed to be a controlling factor.

Differences existed in the type and explanatory degree of factors that affect the spatial distribution of water yield in different geomorphological zones (Figure 4). For example, the explanatory degree of elevation was not significant in platforms, but was 26.1% in hilly areas; the explanatory degree of land use types was only 11.9% in medium relief mountains but was as high as 53.3% in the platform areas. Similar to the results in the whole region, climate-related factors generally had the strongest explanatory degree. Precipitation and actual evapotranspiration contributed the most, and the temperature factor can be ignored. Of the land use related factors, the contribution of land use composition and structure was limited to less than 10% for each factor. In flat terrains such as plains and platform areas, the actual evapotranspiration and land use type were the controlling factors for the distribution of water yield and had a significant difference with other factors. This is likely because the precipitation is relatively uniform on such a small scale causing the actual evapotranspiration and vegetation types to become the principal factors. With an increase in the degree of relief of the land surface, actual evapotranspiration became the controlling factor in hilly areas, and the interpretation power of precipitation has improved and significantly higher than other factors. With the continuous increase of land relief, the contribution of precipitation gradually increased compared to the other factors. In the small to extreme relief mountains, the q values of precipitation exceeded those of actual evapotranspiration to be the first or second controlling factors, and was significantly higher than other environmental factors; the explanatory degree of land use type—one of the controlling factors in the flat areas, was significantly reduced in high-relief areas.The controlling factors for water yield varied across the climatic zones in the Hengduan Mountain region; however, precipitation, actual evapotranspiration, and land use type were still the main influencing factors (Figure 4). In the Qamdo and Bomê-western Sichuan regions in the plateau climate zone, actual evapotranspiration was the factor with the strongest explanatory degree and had a significantly higher impact than land use type or precipitation. This area belongs to the semi-humid region with relative low and uniform precipitation. Regional evapotranspiration had, therefore, become an important factor, and physiological characteristics such as vegetation type and root depth indirectly affected the distribution of water yield by changing the evapotranspiration coefficient. Although all the regions belong to the mid-subtropical climate, the controlling factor types and explanatory degrees varied. In the northern Yunnan region, precipitation is abundant and has obvious spatial heterogeneity. Specifically, the Three Parallel Rivers Region in the western part of the Hengduan Mountain region is the area with the highest precipitation value. The amount of precipitation decreases in the eastward direction ranging from 600 to 1800 mm. Such strong spatial heterogeneity caused precipitation to become the main controlling factor with an explanatory degree of 62%. This value was significantly higher than the explanatory degree of actual evapotranspiration and solar radiation (27%). The interpretation of soil type was limited to only 13%, and the effects of land use type can be ignored. The most influential factor for the spatial heterogeneity of water yield was actual evapotranspiration in the Sichuan region and the Jinsha river-Chuxiong, Yuxi region. The results in these regions were very different from those found in the northern Yunnan region, and land use type plays an important role. This is mainly because the homogenization of precipitation in such a small-scale region causes the actual evapotranspiration to become the dominant factor for the water yield distribution. In the plateau climate zone, actual evapotranspiration and land use types were the dominant factors that explain spatial heterogeneity. In the mid-subtropical climate zone, precipitation was relatively abundant, and its explanatory degree was significantly higher than that of the actual evapotranspiration.

The above content analyzed the influencing degree of a single factor on the water yield spatial distribution. However, in actual, it is the complex interaction of multiple factors determines the spatial pattern. The interaction detector for the paired factors proved that the influence of the interaction of the factors on the water yield spatial distribution was much higher than that of a single factor in either the entire Hengduan Mountain region or the geomorphological and climatic zones, which presented as two factors enhancement and nonlinear enhancement (Table 2). At the mountain scale and in each sub-region, the interactions among climate factors had the strongest explanatory degree, followed by the interaction between climate factors and land use related factors. The q value for the interaction between precipitation and actual evapotranspiration was the highest and explains 97% of the water yield spatial heterogeneity. This indicates that within the same precipitation (actual evapotranspiration) stratum, the spatial difference in actual evapotranspiration (precipitation) would significantly enhance the spatial heterogeneity of water yield even if the precipitation (actual evapotranspiration) were spatially similar. The second and third most important interactions were different in each of the subregions but involved the interaction of precipitation or actual evapotranspiration with other factors. In different geomorphological zones, the combination of precipitation and land use type had a great impact on the water yield, and the explanatory degree reached a high of 70%-90%. For the plateau climatic regions (Qamdo and Bomê-western Sichuan regions), the interaction of actual evapotranspiration with other factors was the major contributor (with a q value of 60%-70%) of which the interaction between actual evapotranspiration and solar radiation was most significant. In the mid-subtropical climate zone, the interaction between actual evapotranspiration and land use type had the highest explanatory degree in the Sichuan region, and the Jinsha river- Chuxiong, Yuxi region. In the northern Yunnan region, the interaction between precipitation and land use type contributed a significant amount, second only to the contribution of the interaction between precipitation and actual evapotranspiration.

Complex linear or nonlinear relationships exist between water yield and various influencing factors. A risk detector module helps to identify the water yield service level in each stratum for various factors and determine if the differences in water yield among the strata of each influencing factor are significant (with a confidence level of 95%), and to identify the areas with high water yield value by analyzing the relationship of water yield with each factor. We did not consider the factors related to land use composition and structure, which belonged to the land use category of influencing factors, into the running of risk detector, because the previous results of the factor and ecological detectors showed that their contribution to the water yield is limited. Areas with high water yield value varied significantly among different geomorphological and climatic zones in the Hengduan Mountain region (Table 3). The relationship between precipitation and water yield presented similar trends of increasing precipitation with increased water yield depth in the entire Hengduan Mountain region and different geomorphological and climatic zones. The water yield changed monotonically with the change in actual evapotranspiration, which indicated that the increase in actual evapotranspiration would reduce the water yield service to a certain extent in both the entire mountain region and the various sub-regions. It seems that the effect of temperature on the spatial heterogeneity of water yield service was more complicated when compared with other climate related factors. In all subregions except for the Sichuan climatic region, water yield showed an increasing trend as the NDVI value increased, which was probably a result of the dense vegetation that was able to intercept more precipitation. In terms of soil type factor, the high values of water yield service were mainly located in the Primarosol and Anthropogenic soil areas. Compared with other land use types, the water yield in construction land was the largest, but it is worth noting that the water yield generated in this area eventually flows into urban drainage pipes and is, consequently, difficult to access. In terms of vegetation area, shrubland had a high water yield value. In the vertical direction, except for the platforms, hills, and the northern Yunnan climatic zone, the water yield gradually decreased with increasing elevation. The decrease may be attributed to the vertical heterogeneity of land use types. The impact of slope on the spatial heterogeneity of water yield was more complicated and inconsistent and might be the result of specific environmental influences in different sub-regions.

Identification of the spatial heterogeneity of ecosystem services and their controlling factors is an important component of ecosystem service research and is the scientific basis for territorial spatial planning and regional resources and environmental carrying capability. The Hengduan Mountain region is located in the upstream area of many major rivers in China and Southeast Asia. The complex topography and diverse climate types have enhanced the complexity of the spatial heterogeneity of water yield service to a certain extent. In this study, we constructed a research framework for the attribution analysis of the spatial heterogeneity of water yield in different geomorphological and climatic zones in mountain regions. Based on the simulated water yield service from the InVEST software, we identified the main controlling factors for the spatial heterogeneity of water yield service and the important areas of water yield by using the geographical detector method. The main conclusions are as follows:

(1) Compared with other types of influencing factors, climatic factors had the strongest explanatory power for the spatial heterogeneity of water yield. Precipitation and actual evapotranspiration were the main factors, and the explanatory power of temperature was extremely limited.

(2) The explanatory power of different factors for water yield spatial heterogeneity varied across different geomorphological zones. In flat areas (including plains and platforms), the actual evapotranspiration had the strongest explanatory power, followed by the land use type. Future regional development should pay more attention to the transfer and changes in the type of land use. The explanatory power of precipitation was gradually enhanced as the degree of relief increased. This area was mainly affected by climatic factors, and the explanatory degree of land use related factors was significantly reduced.

(3) The controlling factors and the explanatory degree of each influencing factor for the spatial heterogeneity of water yield varied depending on climate type. In plateau climatic regions (Qamdo and Bomê-western Sichuan), the actual evapotranspiration had the strongest interpretation ability, and the interactions between this and other factors were the main contributor. In terms of the mid-subtropical climate zone (mainly including the northern Yunnan region), the explanatory power of precipitation increased and exceeded that of actual evapotranspiration to become the main controlling factor. At the same time, the interactions between land use type and other factors had strong explanatory power; therefore, more attention should be paid to regional land use changes.

(4) The distribution of important areas of water yield service was identified by using the risk detector. The Primarosol and Anthropogenic soil areas, shrub distribution areas, and the areas with slopes of less than 5° and between 25° and 35° represented the most important water yield areas. The elevation of high water yield service values varied greatly in different geomorphological and climatic zones. For plains, platforms, small relief mountains, medium relief mountains, and large relief mountains, the high water yield areas were located mainly between elevations of 2069 m and 2484 m. Maximum yield elevation was 3377-3797 m for hilly areas and 2924-3377 m for the extreme relief mountains. In each climatic zone, the distribution of important water yield areas was strongly focused on the lower elevations of the Bomê-western Sichuan region and the Sichuan region, ranging from 347 m to 1591 m. Important water yield elevations in the Jinsha river-Chuxiong, Yuxi region was 1591-2069 m, and the distribution in northern Yunnan and Qamdo regions were relatively higher, ranging from 3797-4548 m and 3797-6808 m.

Ecosystem service is affected by multiple influencing factors and the interactions between different factors. However, the coupling mechanism of multiple driving factors and the contribution of each factor to different ecosystem services remains a challenging problem. In this study, we attempted to identify the main controlling factors for the spatial heterogeneity of water yield service by using the geographical detector method. Both the single factor detector and the analysis of interactions between different factors were consistent with other studies and showed that climate related factors are the main controlling factors and have the strongest explanatory power for the spatial heterogeneity of water yield (Delphin et al., 2016; Sun et al., 2011). The precipitation and actual evapotranspiration were the main controlling factors in different geomorphological and climatic zones (Figure 3). Although the results highlighted the driving roles of climate related factors, the contribution of other factors to water yield service should not be disregarded. For example, in plains, platforms, and hilly areas, land use type was the main influencing factor and was second only to precipitation and actual evapotranspiration overall. The contribution of land use type exceeded 50%, particularly in the plains and platforms (Figure 4). These regions also contain the majority of terrain occupied by cultivated land, cities, and other infrastructure projects related to the implementation of the “National New Urbanization Planning” (Lu and Chen, 2015) and “Rural Vitalization” (Liu, 2018) policies. The distribution of land use types will change accordingly in the future; consequently, future urban and economic development processes should pay special attention to intensive land use and reasonable planning to minimize the impacts on water yield service in mountainous areas. Some studies have shown that landscape patterns such as aggregation or fragmentation have certain promoting or inhibiting effects on the expression of ecosystem services (Su et al., 2012; Jordan et al., 2005). Nonetheless, in this study, the explanatory power of each landscape metric on the spatial heterogeneity of water yield was extremely limited (the q value was all less than 10%). This inconsistency may be related to the scale of the calculation. On the one hand, our study calculated the water yield, landscape metrics, and attribution analysis based on a scale of 1 km, which does not consider the most detailed information on spatial differences in land use structure. On the other hand, the selection of the sliding window size in the calculation of the landscape metrics also affects the research results to a certain extent. Subsequent studies should be microscopic in scale and execute research on the impact of landscape structure on ecosystem services with a special focus on the clarification of the scale effect of land use structure.

A certain degree of uncertainty existed in the entire process of evaluation of water yield service. Firstly, the raster formed meteorological data would be affected by some factors such as the selection of meteorological stations and interpolation methods. The impacts of these uncertainties would be more complex in mountain regions. Secondly, the water yield module in the InVEST software failed to include the impact of complex terrain. Thirdly, we analyzed water yield service for only one year and ignored the intra-annual and inter-annual variations. Future research should focus more on the evaluation and analysis of the time scale for water yield service. Revealing the impacts of land use structure on the spatial heterogeneity of water yield service is a challenge. In future studies, attribution analysis of water yield should be performed on a variety of spatial scales, which would help clarify the scale-dependent issues of the attribution analysis.