The study region covers the northeastern economic zone, the northern coastal economic zone, the middle reaches of the Yellow River and the northwestern economic zone, including 10 provinces (Hebei, Shandong, Heilongjiang, Jilin, Liaoning, Shanxi, Shaanxi, Henan, Gansu, Qinghai), three autonomous regions (Inner Mongolia, Ningxia, Xinjiang), and two province-level municipalities (Beijing, Tianjin). These zones differ at the provincial level in their stages of economic development, manufacturing activities (light industry vs. heavy industry), and the status of available ecological resources (abundant vs. insufficient).

Precipitation data from 756 national meteorological stations in China were acquired from the China Meteorological Data Service Center website (http://data.cma.cn/en). The site-level data were interpolated with the ANUSPLIN software package, which provides an interpolation of noisy multi-variate data using the thin plate smoothing splines from the meteorology station observation data. The quality of the interpolated meteorological dataset has been fully evaluated (error rate less than 7%) (Yu et al., 2004).

Land use data from 1980 to 2010 (in the years of 1980, 1990, 1995, 2000, 2005 and 2010) were derived from Landsat TM/ETM remote sensing images and were generated by manual visual interpretation. The ecosystems were divided into five primary types of farmland, forest, grassland, water body, and desert. The datasets can be found at the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http://www.resdc.cn). ESV in the periods of 1982-1989, 1990-1994, 1995-1999, 2000-2004, 2005-2009, and 2010 were calculated based on land use data in the years of 1980, 1990, 1995, 2000, 2005 and 2010, respectively.

Productivity data in this study were from Community Land Model (CLM4.0), Lund-Potsdam-Jena Dynamic Global Vegetation Model (LPJDGVM), Lund-Potsdam- Jena General Ecosystem Simulator (LPJGUESS), and Organizing Carbon and Hydrology in Dynamic Ecosystems Land Surface Model (ORCHIDEE). The input climate data came from CRU- NCEP (http://dods.extra.cea.fr/data/p529viov/cruncep/). The spatial and temporal resolutions of the process-based productivity model data are 0.5 degrees and 30 days, respectively. The spatial resolution of the ESV results was the same as that of the process-based productivity model data.

CLM4.0 is currently one of the most well-developed land surface process models in the world, which is also the land model for the Common Earth System Model (CESM) and the Common Atmosphere Model (CAM). The surface module has been coupled with multiple climate models, mainly including four parts: biogeophysics, hydrological cycle, biogeochemistry, and dynamic vegetation (Lai et al., 2014). LPJDGVM is a coupled dynamic biogeographic earth model that combines process-based terrestrial vegetation dynamics, terrestrial atmospheric carbon and a water exchange framework, and takes some key ecosystem processes into account, such as photosynthesis, carbon allocation, mortality, resource competition, fire disturbance, and soil heterotrophic respiration (Sitch et al., 2003). LPJGUESS simulates the vegetation dynamics and the mass and energy exchanges between the biosphere and the atmosphere. Based on climatic input data, it provides vegetation-related information such as vegetation types and leaf area indexes for climate systems and land surface programs (Smith et al., 2001; Tang et al., 2014). ORCHIDEE calculates the exchange fluxes of carbon dioxide, water and heat with the atmosphere at an hourly scale, and the daily changes of water and carbon sinks. The model is constructed as two coupled modules. The first module describes the energy and moisture exchange between the atmosphere and the vegetation canopy in a 30-minute time step. The second module simulates the terrestrial carbon cycle, including photosynthesis, respiration, carbon distribution, litter decomposition, and phenology at daily time steps (Krinner et al., 2005; Piao et al., 2006; Tan et al., 2010). These models have been widely used to investigate regional and global terrestrial vegetation growth (Mao et al., 2013; Poulter et al., 2013), and extensively validated against observations across different ecosystems and regions, including China (Piao et al., 2006; Tan et al., 2010; Tao et al., 2010; Piao et al., 2013).

The eddy covariance flux approach based on the theory of micrometeorology enables the direct measurements of functions and processes at the ecosystem scale, such as NEP, energy balance, and greenhouse gas exchange, with a spatial resolution of about 100-3000 m and higher temporal resolution (Baldocchi, 2003). According to the geographical location and ecosystem type, the study region was divided into the typical steppe in northern China, alpine steppe in northwest China, farmland ecosystem in central China, and forest in northeast China. We have obtained eddy covariance flux-based NEP records from ChinaFLUX of four sites representing these four ecosystems from 2003 to 2008. The sites are the alpine grassland ecosystem in Qinghai, the forest ecosystem in Changbai Mountain in Jilin, the temperate typical steppe ecosystem in Inner Mongolia, and the Yucheng farmland ecosystem in Shandong. According to the positions of the eddy covariance flux observation stations, we extracted the simulated NEP data from the process-based productivity model to measure the ratio with the actual NEP data from the eddy covariance fluxes. The closer the ratio is to 1, the better the simulation effect of the model.

Based on the equivalence coefficient table developed by Xie et al. (Xie et al., 2003, 2015), the equivalent factor method (Table 1) was applied to estimate the ESV for farmland, forest, grassland, desert, and water body ecosystems in northern China.

Ecosystem services were divided into 11 categories according to a similar classification system from Xie et al. (2015). A positive correlation generally exists between productivity and service of ecosystem food production, raw material production, air quality regulation, climate regulation, waste treatment, maintenance of soil fertility, habitat services, and cultural functions. Similarly, water supply, regulation of water flows, and erosion prevention have proven to be closely related to precipitation (Li, 2010; Pei, 2013; Xie et al., 2015; Mu, 2016). It should be noted that erosion prevention is also influenced by soil and slope, but we just consider rainfall as the main factor in this paper. Therefore, based on the above-mentioned information and spatial-temporal productivity and precipitation data, we used the following formula to further modify the unit value of ecosystem services dynamically:

where Q_i refers to modified ESV per unit area; P_i is standard ESV per unit area, which is calculated from the equivalent coefficients of ecosystem services and standard equivalence factor that is defined as the economic value of natural grain output per unit area of farmland (Xie et al., 2017), i = 1, 2, …, 11, represents the 11 categories ecosystem services; b refers to the productivity or precipitation on a pixel scale, and B is the average productivity or precipitation of the study region. When i refers to food supply, raw material supply, air quality regulation, climatic regulation, waste treatment, maintenance of soil fertility, habitat services or cultural services, b and B refer to productivity in this calculation. When i refers to water supply, regulation of water flows or erosion prevention, b and B refer to precipitation in this calculation. Data processing and statistical analyses were performed using ArcGIS 10.2 and MATLAB 14.0, respectively. Figures were drawn using ArcGIS 10.2 and OriginPro 2018.

The comparison of NEP values from process-based productivity models and the eddy covariance approach showed that the LPJDGVM model has best performance for the NEP simulation, with the ratios of simulated NEP and observed NEP between 0.94-2.85. CLM4.0 showed higher simulation precision in the Yucheng farmland ecosystem, and LPJDGVM has NEP values closest to the NEP from the eddy covariance approach in the Changbai Mountain forest ecosystem, Inner Mongolia typical steppe and Haibei alpine grassland ecosystem (Table 2).

The total ESV per unit area, which was modified by precipitation and NEP, showed obvious variations among models and regions. At the provincial scale, Jilin, Heilongjiang, Liaoning, Shaanxi, and Beijing had the highest ESV per unit area, while Gansu, Qinghai, Xinjiang and Ningxia had relatively low ESV per unit area.

Figure 1 shows the spatial pattern of ESV per unit area of the different provinces. Generally, ESV decreased from northeast to southwest. The total ESV per unit area modified by the CLM4.0 model and precipitation ranged from 17069 yuan ha^-1 (Ningxia) to 114568 yuan ha^-1 (Heilongjiang), with an average value of 60159 yuan ha^-1 (Fig. 1a). The lowest and highest total ESV per unit area modified by the LPJDGVM model were in Xinjiang (10811 yuan ha^-1) and Heilongjiang (113257 yuan ha^-1), respectively (Fig. 1b). For total ESV per unit area modified by the LPJGUESS model (Fig. 1c), the lowest value of 14487 yuan ha^-1 was in Xinjiang, which is the same as the simulated result from LPJDGVM, while Tianjin had the highest ESV per unit area (91133 yuan ha^-1). Finally, the total ESV per unit area modified by the ORCHIDEE model was highest in Liaoning province (103981 yuan ha^-1), and lowest in Xinjiang (14058 yuan ha^-1), with an average value of 62978 yuan ha^-1 (Fig. 1d).

According to the simulation precision of NEP, the LPJDGVM, LPJDGVM, LPJDGVM and CLM4.0 models were selected as the benchmark models for typical steppe in northern China, alpine steppe in northwest China, forest in northeast China, and farmland ecosystem in central China, respectively. The comparisons show that in different ecosystems, ESV per unit area based on the benchmark model was generally lower than those modified by the other models (Table 3).

The total ESV for 11 categories of ecosystem services showed obvious spatial variations, with higher values from the northeast (i.e. Jilin, Heilongjiang, Qinghai, Inner Mongolia, Shaanxi, and Xinjiang) and the lowest values in central and western provinces (e.g. Beijing, Tianjin, Hebei, Shandong, Shanxi, Ningxia). Specifically, Heilongjiang (CLM4.0: 5.13×10¹² yuan, LPJDGVM: 5.07×10¹² yuan, ORCHIDEE: 4.28×10¹² yuan) and Inner Mongolia (LPJGUESS: 4.5×10¹² yuan) had the highest total ESV, while Ningxia (CLM4.0: 0.08×10¹² yuan, ORCHIDEE: 0.1×10¹² yuan) and Tianjin (LPJGUESS: 0.10×10¹² yuan, LPJDGVM: 0.11×10¹² yuan) had the lowest total ESV (Table 4).

Among different productivity models, LPJGUESS and LPJDGVM generally yielded higher total ESV than CLM4.0 and ORCHIDEE, while LPJDGVM and CLM4.0 presented greater ESV ranges than the other models. In order to comprehensively assess the results of all models, we calculated the averaged ESV of the four models in each province. The averaged ESV of each study region show that there is a similar spatial pattern, with higher ESV distributed in northern China, and lower ESV in some western provinces. Heilongjiang (4.57×10¹² yuan) and Inner Mongolia (4×10¹² yuan) had the highest ESV, followed by Qinghai (2.08×10¹² yuan), Jilin (1.99×10¹² yuan), and Shaanxi (1.58×10¹² yuan). Beijing (0.12×10¹² yuan), Ningxia (0.12× 10¹² yuan), and Tianjin (0.1×10¹² yuan) had the lowest ESV.

For different ecosystem services, the regulation services accounted for the highest component (approximately 70%) of the total ESV in all regions, followed by the supporting services (around 20%), while the provisioning services and cultural services were relatively less important (Fig. 2).

Over the past 30 years, total ESV of the 15 provinces, autonomous regions and municipalities showed an insignificant downward trend. In particular, it fluctuated little before 1999, declined sharply in 2003, and then increased gradually to 2010. The results simulated from LPJGUESS and LPJDGVM are relatively higher than those from CLM4.0 and ORCHIDEE (Fig. 3). To assess the spatial pattern of the ESV interannual trend, we conducted a trend analysis of ESV on the pixel scale. It is obvious that the regions with increasing trends of ESV were distributed in northwestern China (Xinjiang, Qinghai, Gansu, and the northeastern part of Heilongjiang), while pixels with decreasing trends of ESV were concentrated in northeastern China (large parts of Inner Mongolia, Jilin, Liaoning, Ningxia) (Fig. 4).

Furthermore, in order to determine the driving forces of temporal dynamics in ESV, trend analyses of NEP, precipitation, and land use cover change were conducted. Precipitation and NEP derived from different productivity-based models showed little fluctuation (Fig. 5), while land cover changed a lot, such as the areas of forest and grassland which decreased by about 0.1×10⁵ km² and 0.6×10⁵ km², respectively (Fig. 6).

The total ESV, ESV per unit area, and temporal dynamics of ESV modified by the four productivity models presented similar spatial and temporal patterns. For example, ESV gradually decreased from east to west, and the consistency of the inter-annual trend of ESV modified by NEP from the four models was relatively high. However, there are still some differences in the numerical simulation results of different models. NEP simulated by CLM4.0 and LPJDGVM are lower than those of LPJGUESS and ORCHIDEE, which suggests that although different productivity-based models incorporate the dynamics of terrestrial vegetation and atmospheric carbon-water (e.g., photosynthesis, carbon allocation, water, and heat exchange flux), the different models focused on different simulation processes. Furthermore, the process-based productivity model that has the value closest to the measured NEP was deemed as the benchmark model for each region. Generally, the ESV from the benchmark model (LPJDGVM) was higher than those in the other models in typical steppe in northern China, alpine steppe in northwest China, and forest in northeast China, while modified ESV of the farmland ecosystem in central China from the benchmark model (CLM4.0) was lower than those from the other models. The desert steppe is located in the ecologically fragile area of China, with less precipitation and lower productivity. The model parameters are difficult to quantify, leading to some deviations in the results. In addition, the lack of data sources for the CLM4.0 and ORCHIDEE models in the desert steppe in northwest China may also be an important reason for the ESV variations among models.

The ESV evaluated by equivalence factor coefficients, which were modified by the CLM4.0, LPJDGVM, LPJGUESS and ORCHIDEE models, displayed similar spatial patterns: the provinces in northern China (e.g., Jilin, Heilongjiang, Qinghai, Inner Mongolia, Shaanxi, and Xinjiang) had higher ESV, while the other provinces in the central area (e.g., Beijing, Tianjin, Hebei, Shandong, Shanxi, and Ningxia) had relatively low ESV. This is consistent with previous studies (Shi et al., 2012; Xie et al., 2017). Possible driving forces for these spatial patterns are the complex interactions among ecosystems, climate change, and human activities. Provinces in northern China usually had the highest ESV, which is consistent with the spatial distribution of vegetation and regional ecological restoration measures, leading to a higher level of NEP. Some studies have mentioned that the greatest contribution to the total ESV in China comes from the forest ecosystem, which also proves the above view (Xie et al., 2017). At the same time, Qinghai and Xinjiang had low NEP because these provinces receive less influence from the East Asian monsoon climate and hence vegetation growth is strongly restricted. However, the vast geographic land of these two regions resulted in a higher level of ESV. In addition, the more developed provinces and municipalities (e.g., Beijing, Tianjin, Hebei, Shandong, and Shanxi) were usually characterized by lower ESV. This probably resulted from the relatively scare vegetation resources and higher proportions of construction land for heavy industries and human demands. Increasing construction land indicates a higher rate of urbanization to a certain extent, which serves as a significant constraint on the ecosystem services (Liu et al., 2021c). This also supports the view of Xie et al. (2017) that there is a negative correlation between ESV per capita and GDP per capita. In summary, the study area can be divided into three types: the type with dense vegetation and high ESV, which emphasizes the important role of vegetation under better hydrological and thermal conditions; the type which has low vegetation coverage and higher ESV but mainly due to vast geographic area; and the type of an economically developed area with the lowest NEP and lowest ESV due to scarce natural resources.

For different categories of ecosystem services, the regulation services provided the highest ESV, mainly including the water flow regulation services and climate regulation services, which were determined by the total area of the vegetation and the equivalence coefficients of different services in these provinces. Generally, the ESV of several provinces with abundant vegetation resources originated from climate regulation, because climate regulation services are the most important services of forest and grassland. Meanwhile, it should be noted that climate regulation services also contributed more to the ESV of Xinjiang, although the land cover of desert was larger than farmland and forest. This situation probably resulted from the fact that the equivalence coefficient of climate regulation services due to vegetation was higher than that of waste treatment from the desert.

The total ESV of the 15 provinces, autonomous regions and municipalities showed an insignificant downward trend with year. The regions with increasing trends in ESV were concentrated in Heilongjiang, Qinghai, and Gansu, while the economically developed parts in this region, such as Beijing, Hebei, Inner Mongolia, Jilin and Liaoning, experienced decreasing trends of ESV, which is consistent with previous findings (Shi et al., 2012; Lou et al., 2019; Liu et al., 2021a; Liu et al., 2021b; Yang et al., 2022). Land use and land cover (LULC) changes significantly affect regional ESV through obvious variations in ecosystem types and distribution patterns (Lou et al., 2019; Chen and Liu, 2021), which also can explain the insignificant decrease of ESV in our study area. For example, Yan and Zhang (2019) analyzed the impact of land cover change on ESV in the Sanjiang Plain of Heilongjiang in the past 60 years and found that as a result of human activities, a large part of the wetlands had been converted into farmland, leading to a decreasing ESV. We found that the area of forest and grassland in the 15 provinces, autonomous regions and municipalities in northern China decreased insignificantly in the past 20 years, which is consistent with previous studies. For example, grassland degradation was quite severe in the agro-pastoral ecotone of northern China as a result of the considerable substitution of grassland by construction land and barren land; and the area change of construction land (increased) and grassland (decreased) in the Beijing-Tianjin-Hebei region was the largest among the land types (Lou et al., 2019). Therefore, the irrational use of forest and grassland resources and expansion of urban construction land are likely the main driving forces for the insignificant downward trend of ESV.

Since 2000, the Chinese government has carried out ecological management and restoration projects, i.e. The Grain to Green Program (GTGP), North Shelter Forest, and The Natural Forest Protection Project, with an aim of increasing the vegetation coverage and mitigating soil erosion. These restorations have achieved major project results (Lu et al., 2018), such as an increased ESV in Ansai District during 1980-2018 (Han et al., 2021). Compared with the decreasing ESV in our study, such differences probably resulted from the variations in study time scales, in which we only focused on the period during 1982-2010. Furthermore, some previous studies have indicated that there was a turning point of vegetation growth in Northern China in 2010, with a decreasing trend of NDVI before 2010 (Di et al., 2021), which would also have a negative impact on ESV. Therefore, while the afforestation policy and ecological restoration measures in northern China have substantially contributed to CO₂ mitigation and ecological environmental improvement gradually in recently years, they have not improved vegetation ESV to a great extent in the previous period in combination with the influence of climate change (Batunacun et al., 2018).