Spatio-Temporal Change of Multiple Ecosystem Services and Their Driving Factors: A Case Study in Beijing, China

: Driven by rapid urbanization, land use patterns have undergone dramatic changes, which have in turn inﬂuenced ecosystem services (ESs). The government has implemented ecological compensation and conservation actions to mitigate this negative impact, especially in metropolises. However, whether these measures will have the desired effect remains unclear. Therefore, under-standing the spatio-temporal characteristics of ESs and their driving factors are crucial for regional development. In this study, we quantiﬁed carbon storage, water yield and soil conservation services based on land use maps. A Geographical Detector (GD) was used to analyze the driving mechanisms of ES changes in Beijing from 1985 to 2020. The results showed that (1) the obvious landscape pattern changes are urbanization, afforestation and cultivated land degradation in Beijing, (2) the three services showed an increasing trend overall, but the changes were different in each period, (3) in general, land use change is the main factor affecting ESs, and the urbanization and afforestation contributed the most. These results suggest that in highly urbanized metropolises, humans can still balance the demands of regional development and ESs reasonable planning. This study highlights the importance of afforestation for ESs, the necessity of harmonizing environmental concerns and human activities, and the need to conduct ecological management in Beijing to protect the ecological environment and coordinate regional development.


Introduction
Ecosystem services (ESs) are the benefits and goods that people obtain from ecosystems [1], and the ESs concept incorporates the well-being provided by an ecosystem into sustainable management policies [2][3][4]. As an important link between humans and ecosystems, ESs are critical to the sustainable development of human society and the stability of ecosystems [5,6]. However, the pressures of urbanization and global change mean that the ESs provided by ecosystems are facing unprecedented threats, and at least 15 types of global ESs are in a state of degradation [7,8]. Such deterioration damages the restorative capacity of ecosystems and produces irreversible effects, resulting in adverse outcomes for both human development and the ecological environment [9][10][11].
Human activities and climate change have been considered to be the two main factors driving the provision and trade-offs of ESs [12]. Climate change influences ESs by affecting the biophysical processes of ecosystems. The energy flow, material circulation and information transfer of ecosystems fluctuate due to climate change, which may cause irreversible effects on regional ecosystems [13]. Many studies have shown that global warming intensifies the ecological impact of climate change and poses a nonnegligible threat to human survival [14]. For example, under the dual effect of climate change and human activities, many cities around the world are facing water crises [15].
Furthermore, land use change can directly affect the composition of ecosystems by altering their internal processes, which can significantly impact ESs [16,17]. China has experienced rapid urbanization over the past decades, which has been accompanied by rapid changes in land use patterns that have affected regional ESs [18,19]. On the one hand, these changes have led to a series of ecological problems [18,20], such as soil erosion [21], reduced carbon storage and water yield [22], deteriorating water quality [23], and the destruction of biodiversity [24]. These ecological problems threaten the sustainable development of urban areas and endanger the health of urban ecosystems [20,25]. Assessing and mitigating the impacts of urbanization on ecosystems have become priorities for landscape planning. On the other hand, the government has implemented ecological projects to mitigate the impacts of human activities on the ecosystem. For example, Chen et al. [26] found that China has accounted for 25% of the global net increase in leaf area over the past 17 years, which has significantly reduced the CO 2 content in the atmosphere. These results illustrated that the impact of human factors on ecosystems depends on scientific planning, which makes it particularly important to understand the driving mechanisms of ESs changes [27].
In regions with strong human intervention, the effect of socio-economic factors on ecosystem interventions must not be overlooked [28]. Zhang et al. [29] found that the mismatch between supply and demand of ESs in urbanized areas and rural areas has a large gap, and ESs supply and demand were most sensitive to population density in developed urban areas and artificial land proportion in rural areas. Research on the driving mechanisms of regional ES changes is still at the exploratory stage, some studies have recognized and quantified multiple ESs at different scales, and many studies have analyzed the driving factors of ESs [30,31]. However, most of these studies have focused on the effects of land use change, climate change or their interaction on ESs while ignoring socio-economic factors [5]. Integrating these three types of driving factors can help us systematically study the mechanisms driving changes in regional ESs, which can help to make regional policies and implement the ecological projects.
Beijing, as the largest city in China, is facing increasingly serious environmental problems in the process of rapid urbanization. The government of Beijing has implemented some ecological compensation and conservation actions to alleviate the negative effects of urbanization, such as two phases of "One Million-Mu Plain Afforestation Project" from 2012 to 2020 [32]. However, the specific effects of these measures remain unclear. Thus, clarifying the driving mechanisms of ESs changes is critical to designing policies that reduce the negative impacts of human activities and promote sustainable development in the region [33]. This study focused on Beijing, which is one of the most developed urban areas in China, to analyze the spatio-temporal changes in carbon storage, water yield and soil conservation services over the past 35 years. Additionally, the driving mechanisms of the changes in different ESs were explored. Based on the research results, corresponding policy recommendations were developed to improve the ecological environment in Beijing and provide a reliable reference for the sustainable development of megacities in China.

Study Area
Beijing, which is located between 39 • 26 -41 • 03 N and 115 • 25 -117 • 30 E, is the political, economic and cultural center of China. It has a total area of 16,808 km 2 and is dominated by hills and mountains, with the terrain declining from north to south ( Figure 1). The mean annual temperature in Beijing is 11 • C-13 • C, and the yearly average precipitation is 644 mm in the typical warm temperate semi-humid continental monsoon climate zone. Beijing has experienced rapid development since the beginning of the Reform and Opening-up in 1978, with its Gross Domestic Product (GDP) increasing from CNY 10.88 billion to 3032 billion and the total resident population swelling from 8.715 million to

Data Sources
The land use/cover data were acquired from Landsat TM and OLI (downloaded from the United States Geological Survey, https://www.usgs.gov/, last accessed on 3 February 2022) by random forest classification at a 30 m spatial resolution for the years 1985, 1990, 2005, 2000, 2005, 2010, 2015, and 2020. In this study, we analyzed ESs changes based on changes in land use patterns.
The meteorological data used in this study, including daily precipitation, temperature, humidity, air pressure, sunshine, wind speed, and evaporation data, were obtained from the China Meteorological Data Service Center (http://data.cma.cn/, last accessed on 3 February 2022). We interpolated the meteorological data from 94 stations by ANUSPLIN in R and obtained gridded images of the meteorological parameters as the input for the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model. A digital elevation model (DEM) with a 30 m resolution was obtained from the United States Geological Survey (https://www.usgs.gov/, last accessed on 3 February 2022). Soil texture and depth data with a 1 km resolution were acquired from the Harmonized World Soil Database (HWSD) [34]. The dynamic economic, social, and development statistics in this research came from the Beijing Statistical Yearbook (http://tjj.beijing.gov.cn/, last accessed on 3 February 2022). The population density data was download from WorldPop with a 1 km resolution (http:// www.worldpop.org/, last accessed on 3 February 2022) and the GDP data with a 1 km resolution was downloaded from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/, last accessed on 3 February 2022) [35].

Methods
The modified Morris screening method was used to analyze the sensitivity of parameters required for model calculation. Components of the InVEST model and the Revised Universal Soil Loss Equation (RUSLE) model were jointly applied to measure the ESs of Beijing. The InVEST model was used to evaluate the ESs of carbon storage and water yield, and the RUSLE model was used to evaluate soil conservation in the region. Then, we

Data Sources
The land use/cover data were acquired from Landsat TM and OLI (downloaded from the United States Geological Survey, https://www.usgs.gov/, last accessed on 29 January 2022) by random forest classification at a 30 m spatial resolution for the years 1985, 1990, 2005, 2000, 2005, 2010, 2015, and 2020. In this study, we analyzed ESs changes based on changes in land use patterns.
The meteorological data used in this study, including daily precipitation, temperature, humidity, air pressure, sunshine, wind speed, and evaporation data, were obtained from the China Meteorological Data Service Center (http://data.cma.cn/, last accessed on 29 January 2022). We interpolated the meteorological data from 94 stations by ANUSPLIN in R and obtained gridded images of the meteorological parameters as the input for the Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model. A digital elevation model (DEM) with a 30 m resolution was obtained from the United States Geological Survey (https://www.usgs.gov/, last accessed on 29 January 2022). Soil texture and depth data with a 1 km resolution were acquired from the Harmonized World Soil Database (HWSD) [34]. The dynamic economic, social, and development statistics in this research came from the Beijing Statistical Yearbook (http://tjj.beijing.gov.cn/, last accessed on 29 January 2022). The population density data was download from WorldPop with a 1 km resolution (http:// www.worldpop.org/, last accessed on 29 January 2022) and the GDP data with a 1 km resolution was downloaded from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/, last accessed on 29 January 2022) [35].

Methods
The modified Morris screening method was used to analyze the sensitivity of parameters required for model calculation. Components of the InVEST model and the Revised Universal Soil Loss Equation (RUSLE) model were jointly applied to measure the ESs of Beijing. The InVEST model was used to evaluate the ESs of carbon storage and water yield, and the RUSLE model was used to evaluate soil conservation in the region. Then, we calculated the change values for the three ESs from 1985 to 2020, collected and sorted the meteorological, ecological, and socio-economic driving factors, selected eight driving factors through principal component analysis (PCA), and analyzed the driving mechanism of each service change through the Geographical Detector (GD) method. The structure of our study is shown in Figure 2.
Forests 2022, 13, x FOR PEER REVIEW 4 of 18 calculated the change values for the three ESs from 1985 to 2020, collected and sorted the meteorological, ecological, and socio-economic driving factors, selected eight driving factors through principal component analysis (PCA), and analyzed the driving mechanism of each service change through the Geographical Detector (GD) method. The structure of our study is shown in Figure 2.

LULC Maps
In this study, cloudless growth season images were selected, radiometric calibration and atmospheric correction performed in ENVI 5.3, random forest classification applied to classify the corrected image and land use/land cover (LULC) maps were obtained from 1985,1990,1995,2000,2005,2010, 2015, and 2020. Finally, based on the field data from forest inventories of Beijing, the random sampling method was used to verify the accuracy of the results. In the classification results, all images had kappa values above 0.8, which means that the characterized LULC maps were acceptable.

Sensitivity Analysis
The Morris screening method is a global sensitivity analysis method with high efficiency and accuracy, which adopts the One factor At a Time (OAT) type calculation method [36]. In the Morris method, only one parameter is modified within its range and the change result is reflected in the model outcome y = y (x1, x2, …, xn) [37]. Finally, the elementary effect is used to represent the influence of parameter xi change on the output result and can be calculated as follows:

LULC Maps
In this study, cloudless growth season images were selected, radiometric calibration and atmospheric correction performed in ENVI 5.3, random forest classification applied to classify the corrected image and land use/land cover (LULC) maps were obtained from 1985,1990,1995,2000,2005,2010, 2015, and 2020. Finally, based on the field data from forest inventories of Beijing, the random sampling method was used to verify the accuracy of the results. In the classification results, all images had kappa values above 0.8, which means that the characterized LULC maps were acceptable.

Sensitivity Analysis
The Morris screening method is a global sensitivity analysis method with high efficiency and accuracy, which adopts the One factor At a Time (OAT) type calculation method [36]. In the Morris method, only one parameter x i is modified within its range and the change result is reflected in the model outcome y = y (x 1 , x 2 , . . . , x n ) [37]. Finally, the elementary effect e i is used to represent the influence of parameter x i change on the output result and can be calculated as follows: where y * is the new output after the parameter changes, y is the previous output, and ∆ i is the variation in the parameter i. In this study, the modified Morris screening method was used, which perturbs a parameter with a fixed step and calculates the model output [38]. Additionally, the sensitivity index S was calculated as follows: where Y i+1 is the output after the parameter changes for i + 1th, Y i is the output after the parameter changes for ith, and Y 0 is the initial output after parameter calibration. P i+1 is the percentage change in parameter for i + 1th, and P i is the percentage change in parameter for ith. The sensitivity was classified into four levels: highly sensitive parameter (|S| ≥ 1), sensitive parameter (0.2 ≤ |S| < 1), medium sensitive parameter (0.05 ≤ |S| < 0.2), and insensitive parameter (0 ≤ |S| < 0.05).
In this study, we sorted out the parameters needed for model calculation and their value range (Table 1). We disturbed a parameter with −15%, −10%, −5%, 5%, 10% and 15% of the original value, while the other parameters remained unchanged, to explore the sensitivity parameters of three ESs in 1990, 2000, 2010 and 2020.

Carbon Storage
Based on multiple spatio-temporal LULC maps, the InVEST model was used to assess the carbon storage for the different regions and to estimate the carbon stock changes in each unit [39]. The carbon density of the InVEST model includes the aboveground biomass carbon density, belowground biomass carbon density, soil organic carbon density and dead organic carbon density. The relationship among the different carbon densities and carbon stocks is as follows: where C t is the carbon storage in the region at time t, and A jt is the land use pattern at time t. j = 1, 2, . . . represents the different land use patterns in Beijing. Ca j ,Cb j ,Cs j and Cd j are the carbon densities of the aboveground biomass, belowground biomass, soil and dead

Water Yield
The water yield module in the InVEST model provides estimates based on the water balance method [44], which calculates the actual precipitation minus the actual evapotranspiration in each unit. The result of this method is surface runoff, which does not consider the interactions among surface water, subsurface water and baseflow. The annual water yield Y x was calculated as follows: where Y xj represents the water yield of land use pattern j in unit x; AET xj is the annual actual evapotranspiration of land use pattern j in unit x, and P x is the annual actual precipitation of unit x. Evapotranspiration was calculated according to the improved Penman-Monteith formula [45], which uses meteorological data such as temperature, sunshine, and wind speed. Then, we used the evapotranspiration data for the regional site to interpolate the precipitation spatial distribution map by ANUSPLIN in RStudio, used soil texture data to obtain the plant available water fraction [41] and used precipitation data to calculate the Z parameter [46]. Finally, the DEM data and regional statistical data were used to obtain the regional watershed data.

Soil Conservation
Considering the topography in Beijing is dominated by mountains and hills, the RUSLE model was applied to assess the soil conservation. The RUSLE model calculates the difference between the potential and the actual soil loss to reflect the soil conservation [46], and thus the annual soil conservation of the region was calculated as follows: where SR x reflects the soil retention in unit x. R x is the rainfall erosivity factor in unit x, which was estimated by precipitation data in the region; and K x is the soil erodibility factor in unit x, which was estimated by the Erosion-Productivity Impact Calculator (EPIC) [47]. Soil texture (silt, clay, sand) and organic carbon content data from the HWSD were used in this process. L x reflects the slope length factor and S x reflects the slope steepness factor, which we obtained in ArcGIS 10.4.1 based on the DEM data with a 30 m resolution [48]; C x is the crop/vegetation and management factor in unit x, which captures the effect of surface vegetation on soil erosion reduction; and P x represents the support practice factor in unit x, which has a value between 0 and 1. In this study, C x and P x were collected from local studies [49][50][51].

Identification of Driving Factors
According to relevant studies [4,6,30], our research selected 17 driving factors that included three aspects: ecological, socio-economic and meteorological factors. Among the total driving factors, some showed a high redundancy in PCA, and we removed these factors. Finally, we selected eight driving factors to study the changes in ESs, and these factors are shown in Table 2. In particular, LULC changes were selected to represent the total change of human activities in the past 35 years, while GDP data in 2015 and population density data in 2020 were used to represent the spatial distribution of human activities due to the availability of data and non-repeatability of research contents. GD is a model to describe the spatial effects of the explanatory variables on the interpreted variables, including factor detector, risk detector, interactive detector, and ecological detector [27]. The factor detector is commonly used to represent the explanatory ability of driving factors to independent variables, and its result is represented by q value, in the range of [0,1]. The risk detector is used to determine whether there are significant differences in the mean values of attributes between driver sub-regions. GD has been widely used in ecology and economics. Therefore, we used standardized processing of various ES changes combined with the factor detector to explore the influences of different driving factors on ES changes. We used the quantile method to classify the eight driving factors. Then, we divided the independent driving factors into six categories. The categories of LULC represent the conversion of land use types during the study period, while the categories of other factors are arranged in value order from low to high. Finally, we used the risk detector to calculate the mean value of changes in ESs in different categories of each driving factor and analyze the ESs changes in different categories. These processes were performed using the "GD" package in R [52][53][54].

Land Use Changes from 1985 to 2020 in Beijing
The land use patterns in Beijing underwent significant changes from 1985 to 2020 according to Figure 3. Because the terrain of Beijing is high in the Northwest and low in the Southeast, woodland and grassland are mainly distributed in the west and north, while cropland mainly covers the east and south. The area of construction land in the central area showed an obvious increasing trend. The area results of different land use types in each period showed that woodland, cropland and grassland areas accounted for the most abundant, followed by construction land, water, and bare land (Table 3). Construction land, woodland and water showed an increasing trend, while cropland, grassland and bare land showed a decreasing trend during the past 35 years. Among them, the most significant change was the increase in the ratio of construction land from 5.85% in 1985 to 16.36% in 2020, representing a growth area of 1725 km 2 that was primarily converted from cropland. The biggest change was in woodland areas increased by 823 km 2 , while the ratios of grassland, bare land and cropland decreased, especially that of cropland which experienced the greatest decrease of 31%.

Sensitivity of Parameters
There were slight differences in parameter indexes, while the sensitivity rank of parameters was completely consistent in each period. Additionally, the results of sensitivity analysis showed that parameters have different sensitivities to ESs (Table 4). Carbon density was the most important parameter for evaluating carbon storage, precipitation and reference evapotranspiration determined the water yield, and precipitation and soil erodibility had the most significant effect on soil conservation in Beijing. According to the results of sensitivity analysis, we adjusted the parameters to highly sensitive. For example, in order for the precipitation, potential evapotranspiration and rainfall erodibility factor to better reflect the authenticity of the data, we considered the influence of DEM and slope during interpolation. Meanwhile, we referred to the research of Beijing and its environs as much as possible in the process of carbon density data collection [40][41][42][43].

Sensitivity of Parameters
There were slight differences in parameter indexes, while the sensitivity rank of parameters was completely consistent in each period. Additionally, the results of sensitivity analysis showed that parameters have different sensitivities to ESs (Table 4). Carbon density was the most important parameter for evaluating carbon storage, precipitation and reference evapotranspiration determined the water yield, and precipitation and soil erodibility had the most significant effect on soil conservation in Beijing. According to the results of sensitivity analysis, we adjusted the parameters to highly sensitive. For example, in order for the precipitation, potential evapotranspiration and rainfall erodibility factor to better reflect the authenticity of the data, we considered the influence of DEM and slope during interpolation. Meanwhile, we referred to the research of Beijing and its environs as much as possible in the process of carbon density data collection [40][41][42][43].

Spatio-Temporal Changes in ESs from 1985 to 2020
According to Table 5, the three ESs had a decreased and then increased trend during the past 35 years. The total carbon storage increased from 2.38 × 10 8 tons in 1985 to 2.45 × 10 8 tons in 2020, with an increase of 3%. The total water yield service was affected by meteorological factors and fluctuates greatly. Over the total study period, the water yield was the highest in 1990, with a value of 2.01 × 10 9 m 3 , and the total water yield increased from 1.62 × 10 9 m 3 in 1985 to 1.72 × 10 9 m 3 in 2020, with an increase of 6%. The total soil conservation, with trends of fluctuations roughly the same as water yield service, increased from 1.76 × 10 8 tons in 1985 to 1.96 × 10 8 tons in 2020, representing an increase of 11%. It is worth noting that the lowest values of three ESs occurred in the decade from 2005 to 2015. Although ESs have lower growth rates in total, compared with the lowest value the increase rate of ESs was more obvious (water yield service had the highest increase of 156%). The spatial distribution and changes in ESs in different periods from 1985 to 2020 are shown in Figure 4. The high-provision areas of carbon storage and soil conservation services were mainly located in the northern, western, and southern regions, which were the main areas for woodland and grassland in Beijing; these areas had high carbon storage and low soil loss. The low-provision areas of carbon storage and soil conservation were distributed in the central region, which is the main urban region with low vegetation coverage. However, the distribution of water yield was opposite to that of soil conservation and carbon storage, with high-provision areas located in the central region and low-provision areas distributed in the mountainous region. This result is related to the high evapotranspiration of woodland and grassland, as areas with higher evapotranspiration have lower water production. In addition, the three service values provided by water were relatively low, and that provided by cropland was at a moderate level.

Driving Factor Analysis
The GD results showed that the impacts of the driving factors on ESs changes were diverse ( Figure 5). Four socio-economic factors had a clear influence on the ESs changes compared to other driving factors (the q value was generally greater than 0.2), but their impacts are still different. Among them, the q values of GDP and population in the change of three ESs have little difference and the NDVIm has a more significant impact on carbon storage service change among three ESs. The q values of LULC for carbon storage and water yield both exceeded 0.8, indicating that land use change caused by economic development and human activities was the most obvious driving factor of ESs changes in Beijing. The temperature has a stronger impact on ESs changes than precipitation among two meteorological factors. In addition, the slope has the most significant impact on soil conservation service change compared to other driving factors (the q value was greater than 0.4).
According to the classification results of GD in R, we obtained the mean value of ESs changes in each driving factor category ( Figure 6). LULC was mainly classified according to the conversion of land use type, while other factors were divided into six categories ranging from large to small. We highlighted the maximum and minimum ESs changes for each driving factor. The results showed that the changes of driving factors except LULC had no significant impact on the carbon storage service change, which was similar to the result in Figure 5. In all land use type changes, the expansion of woodland had the most significant impact on the increase in carbon storage, and the expansion of construction land had the most significant impact on the increase in water yield and the decrease in soil conservation. The driving results of the water yield service and soil conservation service showed an opposite trend. Furthermore, water yield increased and soil and water conservation decreased with the increase in socio-economic factors and the decrease in ecological factors from 1985 to 2020. The changes of temperature and precipitation of the water yield and soil conservation services were disordered, indicating that the influence of meteorological factors on the two services was not a simple linear relationship.

Driving Factor Analysis
The GD results showed that the impacts of the driving factors on ESs changes were diverse ( Figure 5). Four socio-economic factors had a clear influence on the ESs changes compared to other driving factors (the q value was generally greater than 0.2), but their impacts are still different. Among them, the q values of GDP and population in the change of three ESs have little difference and the NDVIm has a more significant impact on carbon storage service change among three ESs. The q values of LULC for carbon storage and water yield both exceeded 0.8, indicating that land use change caused by economic development and human activities was the most obvious driving factor of ESs changes in Beijing. The temperature has a stronger impact on ESs changes than precipitation among two meteorological factors. In addition, the slope has the most significant impact on soil conservation service change compared to other driving factors (the q value was greater than 0.4).
According to the classification results of GD in R, we obtained the mean value of ESs changes in each driving factor category ( Figure 6). LULC was mainly classified according to the conversion of land use type, while other factors were divided into six categories ranging from large to small. We highlighted the maximum and minimum ESs changes for each driving factor. The results showed that the changes of driving factors except LULC had no significant impact on the carbon storage service change, which was similar to the result in Figure 5. In all land use type changes, the expansion of woodland had the most significant impact on the increase in carbon storage, and the expansion of construction land had the most significant impact on the increase in water yield and the decrease in soil conservation. The driving results of the water yield service and soil conservation service showed an opposite trend. Furthermore, water yield increased and soil and water conservation decreased with the increase in socio-economic factors and the decrease in ecological factors from 1985 to 2020. The changes of temperature and precipitation of the water yield and soil conservation services were disordered, indicating that the influence of meteorological factors on the two services was not a simple linear relationship.

Comparison with Previous Studies
Urban planning is complex and needs to integrate social, environmental, cultural, and other factors. There are many challenges in the integration of ESs into regional planning, especially in metropolises. The land use change was dominated by urbanization and afforestation from 1985 to 2020 in Beijing, which was similar to the result of Sun et al. [19]. The carbon storage, water yield and soil conservation also showed a trend of first decreasing and then increasing during the past 35 years. This result illustrated that local governments were trying to seek a balance between ecological environment and economic development. For example, the government actively carried out two phases of the "One Million-Mu Plain Afforestation Project" from 2012 to 2020 [55], which has increased 600 km 2 of green space in Beijing and is directly reflected in remote sensing images. At the same time, it is affected by ecological engineering as Beijing has sequestered 4.4 million tons of CO 2 annually since 2010. In addition, the results of the factor detector indicated that ecological factors are generally stronger than meteorological factors in influencing ESs changes. This may be related to regional meteorological characteristics that Beijing has, such as a typical rain and heat synchronization climate with a relatively stable climate and less extreme weather disasters. Additionally, the temperature has a stronger effect on changes in ESs than precipitation. This may be because the main tree species are poplar and oak in Beijing, which are more sensitive to temperature differences and are more drought tolerant. Among all ecological factors, DEM has a significant impact on the changes of all ESs (q value over 0.2), which may be related to the landscape pattern and topographic distribution of Beijing; the construction land was mainly distributed in the plain area, while woodland and grassland were mostly distributed in the mountain area where there was an ecological security barrier in Beijing. The results also showed that the LULC was the most important factor of regional ESs changes, which highlighted the dominance of human activities on ESs changes and was similar to the research of Mashizi et al. [56].
However, most studies have ignored the positive effects of human activities on ecosystems [57]. To explore this impact, we analyzed the specific responses of each driving factor to changes in ESs. The results of the risk detector showed that the expansion of woodland has an obvious impact on the increase in carbon storage. Meanwhile, the fluctuations in other factors did not have a significant impact on the change of carbon storage. This may be explained by the fact that the influence of other factors on regional carbon storage changes was unobvious, which was also correlated to the result of the factor detector. Social-economic factors and meteorological factors contributed to the increase in water yield service, while ecological factors had negative effects on water yield. The InVEST model evaluated the regional water yield based on the water balance method, and the high evapotranspiration of woodland and grassland may result in less water production in regions at higher elevations. It should be noted that slope influenced soil conservation service more than the land use change, and soil conservation has been improved significantly with the increase in slope. It showed regions with low human disturbance also existed, even in a highly developed metropolis, and that these areas may be the major ESs supply area. This result highlighted the necessity and importance of reasonable urban planning. In addition, the expansion of construction land leads to a significant increase in water yield and a decrease in soil conservation. This result also indicated that although ongoing ecological projects have improved regional ESs, the negative impacts of urbanization have not been fully offset, especially in areas of urban expansion. One unexpected finding was the extent to which there may be some correlation between ESs trade-off and driving mechanisms. The risk detector results for water yield service and soil conservation service were almost opposite, and the research of Shen et al. [58] showed that there was a clear trade-off relationship between water yield service and soil conservation service in Beijing. It means that the driving mechanism of ESs changes may indirectly reflect the trade-off and synergy between ESs, and this will be the focus of future work.
The results of GD fully revealed the duality of human disturbance to ESs. We, therefore, recommend the following suggestions to improve the quality of living and achieve sustainable development in Beijing: (1) repair and safeguard the ecological conservation function of Beijing through implementation of ecosystem protection projects such as the "Three-North Shelterbelt Program" [59]. (2) Expand the use range of clean energy to reduce carbon emissions. The government of Beijing has vigorously advocated for the use of clean energy such as solar energy in recent years and has achieved significant results, these efforts represent important measures for reducing urban carbon emissions, as proven by Zhang et al. [23]. (3) Strengthen urban greening, limit large-scale urban expansion, optimize urban water yield systems, improve people's awareness of water conservation, and reduce domestic water waste [60]. (4) Plan for regional land functions and restrict the conversion of land use patterns in some areas, such as lake reclamation, to protect the ecological function within Beijing [61].

Strengths and Limitations
Most previous studies focused on exploring the driving factors of ESs in a particular year, but our methodological framework provided a way to explore the driving mechanism and the impact of each factor on ES changes [6,30]. This method requires straightforward data, is simple and easy to implement, and can be applied for driving factor analysis in other regions. The research revealed the impact of the regional landscape patterns changes due to human disturbance on ESs and has certain implications for the region, providing an important reference for future ecological management and landscape planning in Beijing.
Despite the abovementioned strengths, some limitations exist in this study. The models used in the study have been proved to be reasonable and widely used in ESs assessment [6,20], but also have some limitations. The carbon storage module calculated the historical carbon storage based only on the conversion of land use patterns, which is not sufficiently comprehensive [20]; due to the availability of data, we did not refer to the permanent monitoring sites data and NDVI in the evaluation, which meant that the results could not better characterize the impact of afforestation on regional carbon sinks. The InVEST model used the water balance method to calculate the annual water yield, and this module fails to consider the surface water and groundwater [62]. The soil conservation module mainly relies on the RUSLE, and this model represents an improvement over the Universal Soil Loss Equation (USLE) model but is mainly appropriate for plain areas and areas with gentle slopes. The terrain is dominated by mountains in the northwest of Beijing, therefore, the RUSLE model will generate errors when calculating certain rainfall amounts.
In this research, only three key ecosystem services were selected for evaluation due to data availability and regional ecological characteristics. Therefore, we did not fully characterize the status of ESs in Beijing. Our research is only at the city scale because of the characteristics of regional landscape patterns, which means that our research results may not be applicable to other scales. Finally, in the analysis of the driving factors, eight driving factors were selected, but the range of factors may not have been comprehensive. By combining our driving analysis with the research of Sun et al. [20], we found that as the values of the driving factors change, their impacts on different ESs change, and this difference is related to the trade-off and synergy relationships between ESs. In future studies, we will focus on problems caused by a rapid urban expansion in metropolises, subdivide the land use patterns and analyze the regional landscape structure. Meanwhile, more types of ESs should be analyzed such as water purification, food supply and habitat provision; the trade-offs between different ESs should be assessed to acknowledge the complex connections between them; the changes in ESs and their driving factors should be analyzed at multiple scales; a more comprehensive mechanism should be applied to screen and obtain highly correlated driving factors.

Conclusions
In this study, we used land use maps to calculate the carbon storage, water yield and soil conservation services from 1985 to 2020 in Beijing. Furthermore, we analyzed the influence of driving factors on ESs changes. The results show that land use change was dominated by urbanization and afforestation in Beijing over the 35-year period. In the process of rapid urbanization, carbon storage, water yield and soil conservation services showed a first decreased and then increased trend with changes in landscape patterns. The characteristics of the driving mechanisms for changes in the three ESs are different. Socioeconomic factors, especially land use change, had the greatest impact on ESs changes, while meteorological and ecological factors played different roles for each ES change. It is worth noting that this study explored the positive effects and limitations of human disturbances on ESs changes and demonstrated the necessity of positive disturbances in ecological restoration. Our research provides useful examples for the sustainable development of megacities in China. Moreover, future research should pay more attention to exploring the internal relationship between ESs trade-off and driving mechanisms, which can provide more favorable theoretical support for balancing urbanization and regional sustainable development to realize healthy development.