Changes of Ecosystem Services and Landscape Patterns in Mountainous Areas: A Case Study in the Mentougou District in Beijing

: Land use types have been strongly modiﬁed across mountainous areas. This has substantially altered the patterns and processes of ecosystems and the components of ecosystem services (ESs), and could in turn impact the sustainable development. In the mountainous Mentougou district of Beijing, we explored the changes in land use type (cropland, orchard, forested land, scrubland, grassland, bare land, water bodies, wasteland and built-up land), landscape patterns and ESs as well as their interactions during the past 30 years (1985–2014). The ESs included water yield (WY), carbon stocks (CS) and soil retention rate (SR). The results showed that 23.65% of the land use changed and the wasteland decreased by 80.87%. As for ESs, WY decreased by 47.32% since the year 2000, probably due to the increases in temperature and evapotranspiration. Although the decrease of forested land led to the decrease of CS, the increase of vegetation coverage improved SR. CS decreased by 0.99%from 1990 to 2014, and SR increased by 1.38% from 1985 to 2014. Landscape patterns became fragmented and dispersed, and MPS and CS, SHDI and SR were signiﬁcantly negatively correlated. IJI and CS was positively correlated. This indicated that landscape patterns were highly correlated with ESs. In order to maintain the sustainable development of ESs, we should not only plan land use types, but also consider the rationality of landscape patterns.


Introduction
Ecosystem services (ESs) are the benefits or goods that human beings receive, directly or indirectly, from ecosystems [1]. These transactions, emerging from natural environments, generate the needed conditions for human survival [2,3]. ESs are an important bridge between natural ecosystems and human well-being. Global population growth, with the associated socioeconomic development, is a major driver of the land use change (e.g., through agricultural activities, built-up areas and mining, etc.) and ecosystems (e.g., species diversity, climate change, water quality, etc.). This simultaneously changes the energy flow and material cycle of the ecosystems, thereby affecting ESs and human well-being [4,5]. The Millennium Ecosystem Assessment estimates that more than half of the ESs are fading, and this trend will not slow down [3].
Landscape patterns are an indicator of the ecological process [6]. The interactions between landscape patterns and ecological processes drive the overall dynamics of the landscape and present certain landscape functional characteristics [7,8]. This function is related to human needs and

Study Area
Mentougou district is located in the north-western part (115 • 25 E-116 • 10 E, 39 • 48 N-40 • 10 N) of the Beijing in North China and covers a total area of 1455 km 2 ( Figure 1). This area has a mid-latitude continental monsoon climate, and the climate in the eastern plains and western mountainous areas is very different [28]. The average temperature in the eastern plain is 11.7 • C, and the average temperature in the west is 10.2 • C The average annual precipitation is 520.13 mm (1985-2014 a). The annual average sunshine time is about 2470 h. The precipitation increases gradually from the west to the east. Due to the influence of the mid latitude atmospheric circulation and the monsoon, precipitation changes greatly [29]. In Mentougou, the main native tree species are Quercus variabilis, Pinus tabuliformis Carr., Platycladus orientalis (L.) Franco, Larix gmelinii Rupr, etc. [28]. The forests are generally distributed in the mountainous areas above the altitude of 1000 m. In areas lower than 1000 m, plant communities are dominated by scrubs (e.g., Oxytropis aciphylla Ledeb., Lespedeza bicolor Turcz. and Hibiscus syriacus L., etc.) [30]. The soils are divided into three main soil types (mountain meadow, mountain brown soil and cinnamon soil) ( Figure 2).  In Mentougou, the main native tree species are Quercus variabilis, Pinus tabuliformis Carr., Platycladus orientalis (L.) Franco, Larix gmelinii Rupr, etc. [28]. The forests are generally distributed in the mountainous areas above the altitude of 1000 m. In areas lower than 1000 m, plant communities are dominated by scrubs (e.g., Oxytropis aciphylla Ledeb., Lespedeza bicolor Turcz. and Hibiscus syriacus L., etc.) [30]. The soils are divided into three main soil types (mountain meadow, mountain brown soil and cinnamon soil) ( Figure 2). In Mentougou, the main native tree species are Quercus variabilis, Pinus tabuliformis Carr., Platycladus orientalis (L.) Franco, Larix gmelinii Rupr, etc. [28]. The forests are generally distributed in the mountainous areas above the altitude of 1000 m. In areas lower than 1000 m, plant communities are dominated by scrubs (e.g., Oxytropis aciphylla Ledeb., Lespedeza bicolor Turcz. and Hibiscus syriacus L., etc.) [30]. The soils are divided into three main soil types (mountain meadow, mountain brown soil and cinnamon soil) ( Figure 2).

Land Use Classification and Sub-Watershed Division
In this study, we used data from the Landsat Operational Land Imager (OLI) (bands 1, 2, 4, 5, and 7) remote sensing images for nearly 30 years (1985, 1990, 1995, 2000, 2005, 2010 and 2014) [31,32]. We referred to the land use in Beijing, and the Chinese classification system, and classified the land use, into nine types via supervised classification, using raster maps generated at a 30 m spatial resolution [27,33]. The land use included cropland, orchard, forested land, scrubland, grassland, bare land, water bodies, wasteland (including abandoned lands, bare soil and bare vegetation and sparse vegetation.) and built-up land. We chose 450 point (each land use type had 50 points) representative points distributed throughout the study area using stratified random sampling to examine the accuracy of the interpretation using data from Google Earth (Google Earth high-resolution images) (10 × 10 m) (before 2005) and field investigations (after 2005), all above 85%. We also used the ArcSWAT software to divide the Mentougou into 39 sub-watersheds, which provided the basis for the calculation of the ESs and the landscape patterns ( Figure 3).

Land Use Classification and Sub-Watershed Division
In this study, we used data from the Landsat Operational Land Imager (OLI) (bands 1, 2, 4, 5, and 7) remote sensing images for nearly 30 years (1985, 1990, 1995, 2000, 2005, 2010 and 2014) [31,32]. We referred to the land use in Beijing, and the Chinese classification system, and classified the land use, into nine types via supervised classification, using raster maps generated at a 30 m spatial resolution [27,33]. The land use included cropland, orchard, forested land, scrubland, grassland, bare land, water bodies, wasteland (including abandoned lands, bare soil and bare vegetation and sparse vegetation.) and built-up land. We chose 450 point (each land use type had 50 points) representative points distributed throughout the study area using stratified random sampling to examine the accuracy of the interpretation using data from Google Earth (Google Earth high-resolution images) (10 × 10 m) (before 2005) and field investigations (after 2005), all above 85%. We also used the ArcSWAT software to divide the Mentougou into 39 sub-watersheds, which provided the basis for the calculation of the ESs and the landscape patterns ( Figure 3).

Quantification of Landscape Patterns
To assess changes in the structural characteristics of landscape patterns at the sub-watershed level, we selected the mean patch size (MPS), largest patch index (LPI), interspersion and juxtaposition index (IJI), Shannonʹs diversity Index (SHDI) and aggregation index (AI), to characterize landscape patterns [27,34]. All calculations were extracted from the FRAGSTATS 4.2 software manual.

Methods to Quantify the Ecosystem Services
Considering the local conditions, the feasibility of the assessment method, and the availability of data, we chose the subset of the InVEST model and focused on reporting the services for water yield and carbon stocks. We used RUSLE (universal soil loss equation) model to simulate the soil retention rate and employed the InVEST model to estimate the ESs according to the land use types and biophysical data (digital elevation model, climate regulation, potential evapotranspiration, soil types and soil roughness). The simulation of ESs were based on sub-watersheds. The InVEST model clearly indicated that water yield must be calculated in sub-watersheds, which was convincing and accurate, compared with calculations using raster units [35]. Therefore, in order to analyze the relationship between landscape patterns, the value of the carbon storage and soil retention rate was finally selected to be sub-watersheds.

Quantification of Landscape Patterns
To assess changes in the structural characteristics of landscape patterns at the sub-watershed level, we selected the mean patch size (MPS), largest patch index (LPI), interspersion and juxtaposition index (IJI), Shannon's diversity Index (SHDI) and aggregation index (AI), to characterize landscape patterns [27,34]. All calculations were extracted from the FRAGSTATS 4.2 software manual.

Methods to Quantify the Ecosystem Services
Considering the local conditions, the feasibility of the assessment method, and the availability of data, we chose the subset of the InVEST model and focused on reporting the services for water yield and carbon stocks. We used RUSLE (universal soil loss equation) model to simulate the soil retention rate and employed the InVEST model to estimate the ESs according to the land use types and biophysical data (digital elevation model, climate regulation, potential evapotranspiration, soil types and soil roughness). The simulation of ESs were based on sub-watersheds. The InVEST model clearly indicated that water yield must be calculated in sub-watersheds, which was convincing and accurate, compared with calculations using raster units [35]. Therefore, in order to analyze the relationship between landscape patterns, the value of the carbon storage and soil retention rate was finally selected to be sub-watersheds.

Water Yield
The water yield model in the InVEST model is a simplified hydrologic model without considering the surface, underground and basal flow. Based on the theoretical basis of water balance, the Budyko curve and annual precipitation were used for calculation [36]. The water yield of each space unit was equal to the difference between the annual precipitation and actual evapotranspiration. Therefore, the water yield in this model included not only the surface runoff, but also soil water content, canopy interception and litter water holding capacity.
where AET xj is the actual annual evapotranspiration and P x is the average annual precipitation.
Where Y jx is the annual water yield; AET xj /P x approximates the Budyko curve developed by Zhang et al. [37]. R xj is the Budyko aridity index; K xj is the coefficient of vegetation evapotranspiration; ET 0x is the potential evapotranspiration; and ω x is a dimensionless ratio characterizing the natural climate and soil properties.

Soil Retention Rate
The equation (RUSLE) can predict the actual annual amount of soil erosion [38]. The amount of soil loss will be blocked by vegetation, the reduction of soil loss caused by vegetation is defined as soil retention (SEDR). Soil retention was expressed by the subtraction between the potential soil loss under no vegetation coverage (SEDRET) and the actual soil loss (RUSLE). Soil retention is influenced by rainfall intensity, slope and slope length, thus soil retention cannot effectively explain the contribution of the ecosystem services to soil conservation [39]. In order to analyze real soil conservation function of the ecosystem, we used the rate of soil retention and soil loss under no vegetation coverage (soil retention rate) as a measure of the ecosystem contribution to soil conservation.
where RUSLE is the amount of soil erosion, SEDRET is the potential soil loss under no vegetation coverage, SEDT is the soil loss under vegetation coverage; R is the rainfall erosivity, K is the soil erodibility factor, LS is the slope length-gradient factor, C is the crop management factor and P is the support practice factor.

Statistical Analysis
We used data from seven periods and 39 sub-watersheds to analyze the relationship of landscape patterns and ESs (WY, CS, SR), a total of 273 data. Multivariate analysis was performed to explore the influence of the landscape pattern factors (MPS, LPI, IJI, SHDI and AI) on the ESs. After standardization, a redundancy analysis (RDA) was applied (the largest detrended correspondence analysis gradient length values were all shorter than 3.0). After 499 randomly arranged Monte Carlo displacement tests, the significance of all normalized axis eigenvalues was obtained. We performed RDA and acquired four discriminant functions for ESs and landscape, respectively. The first two functions were extracted and used in bi-plots. The correlation analysis and RDA were conducted using IBM SPSS (version 20, Chicago, IL, USA), and the RDA was performed using CANOCO (version 4.5, Ithaca, NY, USA).

Land Use Types, Weather and Socioeconomic Changes from 1985 to 2014 in Mentougou
The study area experienced 344.10 km 2 of land change, and 23.65% of the land changed from 1985 to 2014 ( Figure 4). Forested land and scrubland accounted for more than 75% of the study area. The dynamic change in land use in the last 30 years, from largest to smallest, it was wasteland (80.87%), cropland (66.97%), bare land (41.63%), built-up land (32.40%), scrubland (27.78%), water bodies (11.80%), forested land (7.18%), orchard (6.93%) and grassland (5.33%). The trend of land use changes was the continuous expansion of built-up land and shrunken of cropland. At the same time, the wasteland decreased significantly. Orchards continued to decrease before 2000, and increased after 2000. Scrubland continued to increase. Bare land, water bodies and forested land decreased, and grassland increased slightly.

Statistical Analysis
We used data from seven periods and 39 sub-watersheds to analyze the relationship of landscape patterns and ESs (WY, CS, SR), a total of 273 data. Multivariate analysis was performed to explore the influence of the landscape pattern factors (MPS, LPI, IJI, SHDI and AI) on the ESs. After standardization, a redundancy analysis (RDA) was applied (the largest detrended correspondence analysis gradient length values were all shorter than 3.0). After 499 randomly arranged Monte Carlo displacement tests, the significance of all normalized axis eigenvalues was obtained. We performed RDA and acquired four discriminant functions for ESs and landscape, respectively. The first two functions were extracted and used in bi-plots. The correlation analysis and RDA were conducted using IBM SPSS (version 20, Chicago, IL, USA), and the RDA was performed using CANOCO (version 4.5, Ithaca, NY, USA).

Land Use Types, Weather and Socioeconomic Changes from 1985 to 2014 in Mentougou
The study area experienced 344.10 km 2 of land change, and 23.65% of the land changed from 1985 to 2014 ( Figure 4). Forested land and scrubland accounted for more than 75% of the study area. The dynamic change in land use in the last 30 years, from largest to smallest, it was wasteland (80.87%), cropland (66.97%), bare land (41.63%), built-up land (32.40%), scrubland (27.78%), water bodies (11.80%), forested land (7.18%), orchard (6.93%) and grassland (5.33%). The trend of land use changes was the continuous expansion of built-up land and shrunken of cropland. At the same time, the wasteland decreased significantly. Orchards continued to decrease before 2000, and increased after 2000. Scrubland continued to increase. Bare land, water bodies and forested land decreased, and grassland increased slightly.  (1985,1990,1995,2000,2005,2010 and 2014) from left to right. BUL refers to built-up land, OC to orchard, CL to cropland, WL to wasteland, WB to water bodies, GL to grassland, BL to bare land, FL to forested land, and SL to scrubland.
During the study period, the wasteland decreased by 80.87% (101.66 km 2 ), of which 66.83% was converted to scrubland and 8.84% was converted to orchards (Table 1). Built-up land has increased by 32.40% (21.11 km 2 ). This increased built-up land are mainly expanded by occupying cropland (15.49 km 2 ), orchard (2.58 km 2 ) and wasteland (2.42 km 2 ). During the study period, the scrubland tended to increase (224.94 km 2 ), mainly through the conversion of other land use types. The amount of forested land showed a decreasing trend (120.46 km 2 ), mainly due to its conversion to scrubland. The water area fluctuated little.  (1985,1990,1995,2000,2005,2010 and 2014) from left to right. BUL refers to built-up land, OC to orchard, CL to cropland, WL to wasteland, WB to water bodies, GL to grassland, BL to bare land, FL to forested land, and SL to scrubland.
During the study period, the wasteland decreased by 80.87% (101.66 km 2 ), of which 66.83% was converted to scrubland and 8.84% was converted to orchards (Table 1). Built-up land has increased by 32.40% (21.11 km 2 ). This increased built-up land are mainly expanded by occupying cropland (15.49 km 2 ), orchard (2.58 km 2 ) and wasteland (2.42 km 2 ). During the study period, the scrubland tended to increase (224.94 km 2 ), mainly through the conversion of other land use types. The amount of forested land showed a decreasing trend (120.46 km 2 ), mainly due to its conversion to scrubland. The water area fluctuated little. The weather and socio-economic trends in the Mentougou were analyzed with time series data between 1985 and 2014 obtained from the statistical yearbook throughout the study region ( Figure 5). The regional weather conditions of the Mentougou showed a warming trend. Temperatures rose steadily over the past 30 years, with an average increase rate of 0.0923 • C/year (p < 0.01). During the same period, the estimated precipitation also increased 2.0069 mm·a −1 . Overall, the weather in the study area showed a warming and wetting trend. The GDP showed a significant exponential growth, with an average increase of 1.  The weather and socio-economic trends in the Mentougou were analyzed with time series data between 1985 and 2014 obtained from the statistical yearbook throughout the study region ( Figure  5). The regional weather conditions of the Mentougou showed a warming trend. Temperatures rose steadily over the past 30 years, with an average increase rate of 0.0923 °C/year (p < 0.01). During the same period, the estimated precipitation also increased 2.0069 mm•a −1 . Overall, the weather in the study area showed a warming and wetting trend. The GDP showed a significant exponential growth, with an average increase of 1.

Changes in Landscape Patterns
LPI generally showed a decreasing trend before 2010, decreasing by 3.82%, and a slight increased after 2010. The MPS was 21.12 ha in 1985, decreasing to 19.45 ha in 2005, and increasing to 20.03 ha in 2014. This showed that the landscape patterns gradually fragmented. IJI and SHDI continued to decrease by 16.07% and 10.71% respectively, from 1990 to 2014. It showed that the landscape patterns showed increasingly typical zonal distribution characteristics, but the diversity of landscape patterns was decreasing gradually. The AI generally showed a decreased trend, indicating that the landscape patterns in the study area were gradually dispersed ( Table 2). The high value of MPS was usually distributed in the central region, usually more than 20 ha. The high value of LPI was located in the southeastern part of Mentougou ( Figure 6). After 1995, the LPI in the western region was increased. Most of the regions had a high value area of the IJI, while the low value areas were mainly distributed in some parts of the central and northern regions, usually below 40%. The distribution of the IJI in most regions showed a trend with first an increase from 1985 to 2005 and then a decrease from 2005 to 2014. The high value of SHDI was mainly distributed in the eastern region, usually more than 1.0. The high value of AI was in the central region.

Changes in Landscape Patterns
LPI generally showed a decreasing trend before 2010, decreasing by 3.82%, and a slight increased after 2010. The MPS was 21.12 ha in 1985, decreasing to 19.45 ha in 2005, and increasing to 20.03 ha in 2014. This showed that the landscape patterns gradually fragmented. IJI and SHDI continued to decrease by 16.07% and 10.71% respectively, from 1990 to 2014. It showed that the landscape patterns showed increasingly typical zonal distribution characteristics, but the diversity of landscape patterns was decreasing gradually. The AI generally showed a decreased trend, indicating that the landscape patterns in the study area were gradually dispersed ( Table 2). The high value of MPS was usually distributed in the central region, usually more than 20 ha. The high value of LPI was located in the southeastern part of Mentougou ( Figure 6). After 1995, the LPI in the western region was increased. Most of the regions had a high value area of the IJI, while the low value areas were mainly distributed in some parts of the central and northern regions, usually below 40%. The distribution of the IJI in most regions showed a trend with first an increase from 1985 to 2005 and then a decrease from 2005 to 2014. The high value of SHDI was mainly distributed in the eastern region, usually more than 1.0. The high value of AI was in the central region.

Changes in Water Yield
The water yield (WY) decreased by 44.20% from 1985 to 2014 (Figure 7). The WY reached the highest value in 2010. This is probably due to the weather in this period (high rainfall and low temperature induced evapotranspiration). In addition, the WY in 2000 was the second highest value of the study period (5.

Changes in Water Yield
The water yield (WY) decreased by 44.20% from 1985 to 2014 (Figure 7). The WY reached the highest value in 2010. This is probably due to the weather in this period (high rainfall and low temperature induced evapotranspiration). In addition, the WY in 2000 was the second highest value of the study period (5.

Changes in Water Yield
The water yield (WY) decreased by 44.20% from 1985 to 2014 (Figure 7). The WY reached the highest value in 2010. This is probably due to the weather in this period (high rainfall and low temperature induced evapotranspiration). In addition, the WY in 2000 was the second highest value of the study period (5.

Changes in Carbon Stocks
The changes of CS in the Mentougou are illustrated in Figure 8. During the study period, CS first increased and then decreased, finally increasing slightly after 2010. From 1985 to 1995, CS increased by 24.61 × 10 4 t (1.44%). This showed that the ecological service function of CS in Mentougou was increased, which was mainly related to the decrease of cropland and wasteland and the increase of scrubland. In 1995In -2000In , 2000In -2010, CS decreased by 13.85  10 4 t (0.80%), 5.09  10 4 t (0.30%) and 2.43  10 4 t (0.14%), respectively. This is mainly because the area of forested land has been decreased. Although the area of scrubland and orchard had increased, but the forested land had the greatest CS compared with other land use types, and its decrease had led to the decline in CS. From 2010 to 2014, the CS increased slightly by 4.53  10 4 t (0.27%). This was closely related to the increase in forested land and the decrease of wasteland in this period.

Changes in Carbon Stocks
The changes of CS in the Mentougou are illustrated in Figure 8. During the study period, CS first increased and then decreased, finally increasing slightly after 2010. From 1985 to 1995, CS increased by 24.61 × 10 4 t (1.44%). This showed that the ecological service function of CS in Mentougou was increased, which was mainly related to the decrease of cropland and wasteland and the increase of scrubland. In 1995In -2000In , 2000In -2010, CS decreased by 13.85 × 10 4 t (0.80%), 5.09 × 10 4 t (0.30%) and 2.43 × 10 4 t (0.14%), respectively. This is mainly because the area of forested land has been decreased. Although the area of scrubland and orchard had increased, but the forested land had the greatest CS compared with other land use types, and its decrease had led to the decline in CS. From 2010 to 2014, the CS increased slightly by 4.53 × 10 4 t (0.27%). This was closely related to the increase in forested land and the decrease of wasteland in this period.

Changes in Carbon Stocks
The changes of CS in the Mentougou are illustrated in Figure 8. During the study period, CS first increased and then decreased, finally increasing slightly after 2010. From 1985 to 1995, CS increased by 24.61 × 10 4 t (1.44%). This showed that the ecological service function of CS in Mentougou was increased, which was mainly related to the decrease of cropland and wasteland and the increase of scrubland. In 1995-2000, 2000-2005 and 2005-2010, CS decreased by 13.85  10 4 t (0.80%), 5.09  10 4 t (0.30%) and 2.43  10 4 t (0.14%), respectively. This is mainly because the area of forested land has been decreased. Although the area of scrubland and orchard had increased, but the forested land had the greatest CS compared with other land use types, and its decrease had led to the decline in CS. From 2010 to 2014, the CS increased slightly by 4.53  10 4 t (0.27%). This was closely related to the increase in forested land and the decrease of wasteland in this period.

Changes in Soil Retention Rate
The spatial distribution of SR was relatively stable, with a higher value usually in the north and south of the mountainous areas ( Figure 9). The SR in the seven periods was 97.35%, 97.83%, 96.00%, 98.40%, 97.93%, 98.72% and 98.70%, respectively. The lowest value for the whole study period was in 1995. The SR was relatively stable after 2000. The highest value was in 2010. Compared with 1995, the SR increased by 2.93% in 2000 and 2.72% in 2014. Generally speaking, the soil conservation function has improved in recent years. Soil conservation increased by 1.35% from 1985 to 2014.

Changes in Soil Retention Rate
The spatial distribution of SR was relatively stable, with a higher value usually in the north and south of the mountainous areas ( Figure 9). The SR in the seven periods was 97.35%, 97.83%, 96.00%, 98.40%, 97.93%, 98.72% and 98.70%, respectively. The lowest value for the whole study period was in 1995. The SR was relatively stable after 2000. The highest value was in 2010. Compared with 1995, the SR increased by 2.93% in 2000 and 2.72% in 2014. Generally speaking, the soil conservation function has improved in recent years. Soil conservation increased by 1.35% from 1985 to 2014.

The Relationship of Ecosystem Services and Landscape Patterns
This paper drafted a three-dimensional coordinate system, in which the XYZ axis respectively represented the changes of carbon stocks, soil retention rate and water yield from 1985 to 2014, as shown in Figure 10. The coordinate points of the levels in ecosystem services (ESs) were distributed on the diagonal body diagonal of the three-dimensional coordinate system. In the past 30 years, the ESs had changed considerably. In general, the SR after 2000 was higher than that before 2000, but the amount of WY and CS showed the opposite result.
We used the landscape indices and ESs of 39 sub-watersheds in seven periods (1985-2014) to carry out the RDA analysis. The overall test of the canonical relationship was significant (p < 0.01). The cumulative percentage variance of ESs explained on the first axis of the RDA was 10.20% (eigenvalues: 0.102), the second axis of the RDA was 1.30% (eigenvalues: 0.013). As shown in Figure  11, the landscape pattern factors related to the height of the first ordination axis were IJI and LPI, while MPS, SHDI and AI were negatively correlated. In the correlation analysis, there was a significant negative correlation between MPS and CS (F = −0.311, p < 0.001), and a significant positive correlation between IJI and CS (F = 0.280, p < 0.001), SHDI and SR were significantly negatively correlated (F = −0.124, p < 0.05). In the ESs of the sub-watersheds, CS was a relatively obvious indicator. MPS, IJI and SHDI were the landscape pattern factors that had great impact on ESs.

The Relationship of Ecosystem Services and Landscape Patterns
This paper drafted a three-dimensional coordinate system, in which the XYZ axis respectively represented the changes of carbon stocks, soil retention rate and water yield from 1985 to 2014, as shown in Figure 10. The coordinate points of the levels in ecosystem services (ESs) were distributed on the diagonal body diagonal of the three-dimensional coordinate system. In the past 30 years, the ESs had changed considerably. In general, the SR after 2000 was higher than that before 2000, but the amount of WY and CS showed the opposite result.
We used the landscape indices and ESs of 39 sub-watersheds in seven periods  to carry out the RDA analysis. The overall test of the canonical relationship was significant (p < 0.01). The cumulative percentage variance of ESs explained on the first axis of the RDA was 10.20% (eigenvalues: 0.102), the second axis of the RDA was 1.30% (eigenvalues: 0.013). As shown in Figure 11, the landscape pattern factors related to the height of the first ordination axis were IJI and LPI, while MPS, SHDI and AI were negatively correlated. In the correlation analysis, there was a significant negative correlation between MPS and CS (F = −0.311, p < 0.001), and a significant positive correlation between IJI and CS (F = 0.280, p < 0.001), SHDI and SR were significantly negatively correlated (F = −0.124, p < 0.05). In the ESs of the sub-watersheds, CS was a relatively obvious indicator. MPS, IJI and SHDI were the landscape pattern factors that had great impact on ESs.

Impact of Land Use Type and Weather Changes on Ecosystem Services
Over the past 30 years, 23.65% of the land use types had changed (Table 1). Of these, the area of scrubland, wasteland and forested land were the main land use types that affected ESs. A large amount of wasteland converted into scrubland, contributing to CS. At the same time, a large amount of forested land changed. Forested land, as the most effective land use type in terms of carbon stock capacity, was one of the reasons for the decrease of CS in the study area after 1995.
In general, the land use types were covered from low vegetation to high vegetation, resulting in increased evapotranspiration and reduced WY [42,43]. This result was also verified in our research.

Impact of Land Use Type and Weather Changes on Ecosystem Services
Over the past 30 years, 23.65% of the land use types had changed (Table 1). Of these, the area of scrubland, wasteland and forested land were the main land use types that affected ESs. A large amount of wasteland converted into scrubland, contributing to CS. At the same time, a large amount of forested land changed. Forested land, as the most effective land use type in terms of carbon stock capacity, was one of the reasons for the decrease of CS in the study area after 1995.
In general, the land use types were covered from low vegetation to high vegetation, resulting in increased evapotranspiration and reduced WY [42,43]. This result was also verified in our research.

Impact of Land Use Type and Weather Changes on Ecosystem Services
Over the past 30 years, 23.65% of the land use types had changed (Table 1). Of these, the area of scrubland, wasteland and forested land were the main land use types that affected ESs. A large amount of wasteland converted into scrubland, contributing to CS. At the same time, a large amount of forested land changed. Forested land, as the most effective land use type in terms of carbon stock capacity, was one of the reasons for the decrease of CS in the study area after 1995.
In general, the land use types were covered from low vegetation to high vegetation, resulting in increased evapotranspiration and reduced WY [42,43]. This result was also verified in our research.
Except for the special weather in Beijing in 2010, there was an abnormally high value of WY. Taking 2000 as a turning point, it presented a trend of first increase and then decrease. This was probably due to the Grain for Green program in the study area in 2000. The increase of scrubland, forested land and orchard was probably due to the increases in evapotranspiration, resulting in low WY.
The SR in the study area showed a gradual increasing trend. This was probably due to the ecological construction and vegetation coverage increaseing in the study area. The increased grassland, scrubland and forested land could increase SR. However, at the same time, this would also increase the evapotranspiration. Moreover, although the scrubland had increased, the forested land was generally decreased. Forested land as the largest carbon reserve led to the reduction of CS. The trade-off of ESs is a question worth pondering. In addition, climate change has a great impact on ESs [36,44]. Abnormal temperature has significant effects on ESs. In the study of the effects of rainfall and temperature on crops, forests and biological activities, researchers have found that warm weather would reduces grain yield, and increases the frequency of forest pests, and the early pollination of bees [45][46][47].

Mechanisms Underlying the Impact of Landscape Pattern Changes on Ecosystem Services
The evolution of landscape patterns is closely related to the change of ESs [14]. However, it is often neglected in many studies of land use types and weather change [36,43,48]. In the study, three indices of landscape indices (MPS, IJI and SHDI) were significantly correlated with CS and SR ( Figure 11). This indicated that the landscape patterns were closely related to ESs. The same conclusion was also confirmed for farming and fertilization, which change the physical and chemical properties of soil, affecting the storage of organic matter and increasing the occurrence of non-point source pollution. Ultimately, this affected the value of the ESs [47,49]. Landscape patterns changes could also be related to the alteration of hydrological cycles, microclimates in the watershed, and soil nutrients [50,51].
Landscape patterns and ESs were significantly correlated [34,52,53]. The positive and negative relationships between these landscape indices and different types of ESs were different in different places. For example, the relationship between landscape patterns and ESs in the South Bay area of Hangzhou Bay indicated that increasing SHDI was conducive to improving the overall ESs [54]. There was a negative correlation between SHDI and ESs in the typical Karst area of Northwest Guangxi. In this study, SHDI was negatively correlated with the soil conservation function. The main reason was that the decrease of SHDI was mainly due to the reduction of wasteland in the study period [55]. The wasteland was mainly transformed into scrubland and forested land with stronger soil conservation capacity, thus SR was improved. In this study, MPS was negatively correlated with CS, which was different from most studies. The main reason for this was a large amount of broken small block scrubland around forested areas. This may result in smaller MPS in scrubland and forested land. It could be seen that the landscape patterns were closely related to the change of ESs and the positive or negative correlation was made according to local conditions.

Strategies for the Sustainable Use of Ecosystem Services
Rational land use planning should take account of the value and sustainable development of multiple ESs and human well-being. This is still a huge challenge and could not really be put into action in many places. In the past 30 years, the built-up land in the study area had increased by 32.40%. We found that 73.38% of these increased built-up lands came from cropland. Moreover, these increased built-up lands were mainly concentrated in the few plain areas in the study area. Generally speaking, people usually choose fertile land for farming. Therefore, the productivity of land had been largely draining away [26,56]. Many ecological control projects have been implemented in the study area, but the overall tendency for forested land was still in decline [57].
Matthew et al. emphasized the impact of fragmentation in landscape patterns on ESs, suggesting that the fragmentation of landscape patterns would fragment biological pathways [52]. For example, the gene exchange of the population is limited by the population of the habitat patches [58,59].
The connectivity and diversity of landscape would also have a significant impact on ESs [34,53]. Our conclusion also proved that the landscape pattern was the main factor affecting ESs.
When people consume one or several ESs, they will intentionally or unintentionally influence the provision of other ESs [60]. This changes the trade-off phenomenon of ESs were changed [61]. It is very important to balance the relationship of this tradeoff. In this study, the increase of scrubland and forested land would have a positive impact on CS. However, it also increased evapotranspiration, resulting in reduced WY. The same trade-off is also reflected in human activities, such as agricultural reclamation and deforestation. These activities, on the one hand, increase services such as food supply and timber supply. On the other hand, it also led to a decline in carbon fixation, soil and water conservation [12]. This contradictory relationship makes the trade-off of ESs increasingly prominent.

Conclusions
Our results showed that 23.65% of the land use types in Mentougou district, a mountainous area in Beijing, changed from 1985 to 2014, with a largest fraction of change being for forested land and scrubland. Overall, the water yield decreased roughly 47.32% after year 2000. Carbon storage generally decreased. The soil retention rate increased by 1.38%. This was probably due to the increase in temperature and evapotranspiration, wasteland reduction and the increase of vegetation coverage. The landscape patterns were dispersed and fragmented. MPS and CS, IJI and CS were significantly correlated (p < 0.001), and SHDI and SR were significantly correlated (p < 0.05). Two important measures should receive attention to maintain the sustainable development of ESs in mountainous areas, (1) rational management of land use types. The loss of basic cropland and the destruction of forested land should be avoided. The influence of evapotranspiration on WY should also be considered during afforestation, and (2) The rationality of landscape patterns should be considered. Landscape patterns affect ESs, but it is often neglected in current models and research. Some landscape pattern indices were significant correlated with ESs.