Spatiotemporal heterogeneity correction in land ecosystem services and its value assessment: a case study of the Loess Plateau of China

The considerable variation in structures and functions of different ecosystems leads to highly variable ecosystem service values (ESVs). Consequently, the accurate quantification of ESVs can better assess and reflect impacts of land use and cover changes (LUCC) on ecosystem services. In the land use simulations of this study, a CA–Markov model was chosen and nine factors affecting land use change were evaluated, followed by the construction of a multi-criteria evaluation method to simulate land use scenarios between 2025 and 2030 on the Loess Plateau. Six key ecological indicators (economy, water production, net primary plant productivity, habitat quality, accessibility, and soil conservation) were used to correct for spatiotemporal heterogeneity within the terrestrial ESV equivalent weight table for China to obtain an ESV equivalent weight factor table that is applicable to the Loess Plateau. Using the newly corrected table, ESVs for the Loess Plateau region were estimated between 1995 and 2030, and the impacts of LUCC on ESVs were analyzed. The Kappa values for the 2015 land use simulation results were 0.80 and 0.83, which were greater than 0.75, indicating that the CA–Markov model simulations were accurate. Throughout the study period, the largest increases in land use type area were for built-up areas and forest lands, with built-up areas primarily derived from conversions of cultivated lands and grasslands, and forest land increases primarily coming from conversion of grasslands. ESVs increased overall by 933.97 × 108 yuan and 312.86 × 108 yuan from 1995 to 2018 and 2018 to 2030, respectively. The three largest contributors to ecosystem services among land use types were moderate grasses, shrublands, and dense grasslands. In conclusion, ESVs for the Loess Plateau steadily increased year by year from 1995 to 2030, indicating that ecological restoration projects played major roles in improving the stability and sustainability of the region’s ecosystems.


Introduction
Ecosystem services refer to the various benefits that humans receive directly or indirectly from environments and are important features of regional ecological civilization and sustainable development (Costanza et al. 1997;Wang et al. 2018;Shi et al. 2021). Ecosystems provide numerous services to humans including supporting, regulating, provisioning, and cultural services (Costanza et al. 2014). The valuation of ecosystem services (ESVs) has become an important research focus (Snyder 2019) that has significantly increased human awareness of ecosystem provisioning well-being and provided a scientific framework to better manage and sustainably develop environmental assets Wang et al. 2016a, b). Human activities have significantly altered much of the Earth's surface in recent years (de Jong et al. 2021). Increasing demands for food, water, and shelter from growing human populations have resulted in largescale deforestation, expansion of urban construction areas (Aburas et al. 2016), fragmentation and loss of agricultural land (Ellis et al. 2013;Long et al. 2014), and accelerated Responsible Editor: Philippe Garrigues * Hengjia Zhang zhanghj@gsau.edu.cn 1 evolution of regional land use patterns globally, all of which have inevitably affected changes in and the sustainable development of land ecosystem services (Chen et al. 2020;Liu et al. 2005). Consequently, the analysis and assessment of ecosystem service distributions caused by LUCC, in addition to their trends, is urgently needed to inform the responsible exploitation of natural resources. Land use and cover change (LUCC) are inextricably linked to ecosystem services, and different land use types exhibit different land coverage that plays important and various roles in ecosystem service provisioning (Mendoza-González et al. 2012;Abulizi et al. 2016;Hu et al. 2020). In addition, LUCC reflects ecological processes and directly causes changes in ecosystem types, areas, and spatial distributions, in turn impacting the structures and functions of ecosystem services . Investigations focusing on the dynamic responses of ecosystem services to land use, ESV assessments, and LUCC modeling have significantly increased in recent years (Azadi et al. 2017;Chen et al. 2020;Liu et al. 2021). Nevertheless, the land use data employed in these studies suffer from low precision and simple (single) classification (Shi et al. 2021;Song and Deng 2017;Wang et al. 2018). In contrast, high precision and multi-class land use data are more appropriate for these studies because they will produce more rigorous assessments. In addition, a significant factor that influences the precise assessment of ESVs is whether the framework used for valuation is reasonable.
Functional value and equivalent factor assessments are currently employed ESV assessment methods. However, the functional value assessment method's high implementation costs; difficulty in obtaining, calculating, and setting ecological parameters; and uncertainty in ESV results render it suboptimal. Costanza et al. (1997) first developed an ecosystem assessment framework for the equivalent factor assessment method that was subsequently used to quantify ESVs at different scales (Costanza et al. 2014;Sannigrahi et al. 2018). However, the application of this ESV assessment model to China has sparked considerable debate. For example, the ESV unit reflecting cropland was grossly undervalued and did not consider a Chinese-specific ESV. In particular, it is difficult to directly monetarily quantify the value of regulatory and support services, leading to several ESVs not being considered or evaluated. Xie et al. (2015) combined expertise and meta-analysis to scientifically and systematically determine an equivalence factor table for China based on Chinese ecosystems and socio-economic development. The factor table has been subsequently widely applied in many Chinese regions (Abulizi et al. 2016;Liu et al. 2021;Xie et al. 2015;Xue et al. 2018). However, the table only statically reflects average ecosystem service function values in China and does not consider spatiotemporal differences of ecosystem types, quality conditions (Han et al. 2020), and dynamic changes. For example, functions related to raw material production, gas regulation, and nutrient cycling (Cao et al. 2017) are closely related to net plant primary productivity (NPP). Further, functions like water supply and hydrological regulation, soil formation and retention, biodiversity protection, recreation, and culture are related to precipitation and evapotranspiration (Sharp et al. 2018), soil erosion (Sun et al. 2018), habitat quality reachability (Paracchini et al. 2014), respectively. Moreover, the ESV equivalent factor table developed by Xie et al. (2015) used the average net profit per hectare value for three major crops (wheat, rice, and maize) in 2010, which clearly does not reflect current realities. Therefore, an accurate and dynamic ecosystem service valuation framework is urgently needed to achieve corresponding spatiotemporal corrections to the equivalent factor table and accurately measure regional ecosystem service capacity.
To address the abovementioned problems in current land ecosystem service valuation, this study employed land use data from 1995-2018 (reclassified into 11 categories) at a resolution of 30 m, using the Loess Plateau of China as an example. A multi-criteria evaluation method was constructed in this study and used with the cellular automata (CA)-Markov model to simulate land use scenarios in 2025 and 2030 for the Loess Plateau. Additionally, six key ecological indicators including economy, water yield, net primary productivity of plants, habitat quality, accessibility, and soil conservation were selected to correct for the spatiotemporal heterogeneity of the Chinese terrestrial ESV equivalent weight table, and a dynamic and accurate ecosystem service value assessment method was subsequently obtained. We hypothesized that the evaluation method could be well used to accurately assess the ESVs of the Loess Plateau from 1995 to 2030. Consequently, the objectives of this study were to determine: (1) the land use structure, changes, and land use interconversions in the Loess Plateau region between 1995 to 2030; (2) a revised ESV assessment framework to assess ESVs in this region between 1995 and 2030 while analyzing the impact of LUCC on ESVs; and (3) the spatiotemporal responses of LUCC to ESVs and temporal responses to ecosystem service functions.

Study area
The Loess Plateau (33° 41′-41° 16′ N, 100° 52′-114° 33′ E) is one of the four major plateaus of China, sitting at an altitude of 12-5,231 m and comprising a total area of about 64.9 × 10 4 km 2 . The plateau features the most concentrated and largest loess areas of the world (Fig. 1), while also being one of the most critically fragile areas of the world due to soil erosion and environmental degradation (Huang and Shao 2019;Zhao et al. 2013). The Loess Plateau contains loose soils and much of its terrain is fragmented with hilly gullies. Rainfall on the plateau mostly occurs in the form of torrential downpours from July to September. The overall climate of the region is dry, with average annual precipitation and evapotranspiration of 490 mm and 350 mm, respectively, in addition to average annual temperatures of 3.6 to 14.3 °C. The climate type of the plateau is classified as an Asian semi-arid continental monsoon climate zone. The Loess Plateau area is the main hydrologically flowing region of the Yellow River and is the largest source of sediment for the lower reaches of the Yellow River, providing nearly 90% of the sediments for the lower reaches (Wang et al. 2016a, b).

Data sources
The land use data (1995, 2000, 2005, 2010, 2015, and 2018), net plant primary productivity (NPP) data (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010), and precipitation data (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015) used in this study were obtained from the Data Center for Resource and Environment Sciences at the Chinese Academy of Sciences (https:// www. resdc. cn/), with a spatial land use data resolution of 30 m, with the other spatial resolutions being 1 km. Land use data derived from Landsat TM/ ETM images were generated through visual interpretation based on national field surveys. Such data have been previously used in similar studies and exhibit an accuracy of over 90% (Chen et al. 2020). NPP data were obtained by calculating the light energy utilization model GLO_PEM, and the annual precipitation spatial interpolation dataset was generated with daily observation data from over 2,400 meteorological stations nationwide through collation, calculation, and spatial interpolation processing. The meteorological data were generated as components of national field surveys. Elevation data were obtained from the Geospatial data cloud (https:// www. gsclo ud. cn/), and slope data were calculated based on elevation data using the ArcGIS software program. Potential evapotranspiration data (2010-2014) were downloaded from the Global Aridity and PET Database (www. cgiar-csi. org/ data/ global-aridi ty-and-pet-datab ase) and exhibit a spatial resolution of 1 km. The area and cost-return data for the three crops that were the focus of the study were obtained from Road network data were obtained from OpenStreeMap (https:// www. openh istor icalm ap. org/). Distances to railroads, provincial roads, and national roads were obtained based on Euclidean distance calculations using ArcGis with a spatial resolution of 300 m. Distance from town construction, rural settlements, and deserts were obtained using ArcGis reclassification of land use data for each time period that was considered and Euclidean distance calculations at a spatial resolution of 300 m. The spatial coordinate systems used in this paper were all Krasovsky_1940_Albers.

Land use simulation
Land use simulation has become an intense area of research, with the main land use simulation methods including system dynamics, Markov models, Agent, CLUE-S simulations, CA, and CA-Markov models. CA-Markov models integrate the features of CA and Markov theory for temporal and spatial predictions, thereby enabling better spatiotemporal simulations of land change from quantitative and spatial perspectives (Sang et al. 2011). For the above reasons, CA-Markov modeling was used in this study to simulate future land use scenarios along with nine mechanistic factors underlying land use characteristics including elevation, slope, rainfall, distance from national roads, distance from provincial roads, distance from railroads, distance from deserts, distance from rural areas, and distance from cities. The factors were also adjusted based on the actual characteristics of the study area ( Fig. 2), with the specific operation steps and CA parameter settings shown in the Appendix.

CA-Markov model
The expression formulas for the CA-Markov model (Sang et al. 2011) are as follows: , and (i) indicate elevation, slope, and precipitation, in addition to distances from national roads, railways, provincial roads, urban, rural areas, and deserts in the study area, respectively where S (t) and S (t+1) are the states of the landscape structure at t and t + 1 times, respectively, while P ij denotes the transfer state.
The expression formulas for P ij are as follows: where t and t + 1 denote the times before and after the cellular automata, respectively; S is the discrete, finite set of states of the cellular automata; N is the domain of the cellular automata; and f is the spatial transformation rule for the cellular automata.

Kappa test
The expression formula for the Kappa Test model is as follows: where P 0 is the proportion of correct simulations; P c is the expected proportion of correct simulations in random scenarios; P p is the amount of correct simulations in the ideal classification case (i.e., 100%).

Degree of land use dynamism (K)
The expression formula of K  is as follows: where K is the degree of land use dynamism during the study period; U a and U b are the areas of a land use type at the beginning and end of the study period, respectively; T is the study period in years.

Calculation of ESVs
ESV was calculated using the following equations: (1) where ESV m denotes the ecosystem service value for landuse type m ( m = 1, 2 ⋯ 12 ); ESV f denotes the value for ecosystem service function type f ( 1, 2 ⋯ 11 ); V i,f denotes the modified ESV equivalent factor table for area i ; A m is the area (ha) for land-use type m ; D r indicates the value of ecosystem services for a modified 1 standard equivalent factor (yuan/ha).
The formula for the service functions is as follows: where V f denotes the average Chinese ecosystem service value equivalent factor scale, and R (i,k) denotes the corrected factor of type k ( k = 1, 2 ⋯ 6 ) in region i.

Economic spatiotemporal correction factor
The national ESV of 1 standard equivalent factor was corrected in this study with the average gross production value for the three major crops, followed by correction of a regional ESV equivalent factor table with the proportion of the three major crops grown relative to the gross production value. The economic value correction formulae were defined as follows: where D 2010 and C 2010 denote the average net profit and total value of output per hectare for the three crops in 2010, respectively; C n denotes the average total value per hectare for the three crops in 2013-2018; E indicates the value of ecosystem services for a modified 1 standard equivalent factor (yuan/ha) in the study area; S r , S w , and S c denote the area sown for rice, wheat, and maize as a percentage (%) of the total area sown for the three crops in the study area in 2013-2018, respectively; and A r , A w , and A c denote the total nationwide value of yield per hectare (yuan/ha) for rice, wheat, and maize from 2013 to 2018, respectively.

Water production spatiotemporal correction factor
Water production affects ESVs by influencing water supply and hydrological regulatory services. The contribution of water supplies from different landscapes can be estimated using the InVEST model that reflects how changes in land use patterns affect annual water production (Sharp et al. 2012). Land use structure also varies from region to region and impacts hydrological processes such as evapotranspiration. Consequently, spatial and temporal corrections for water supply and hydrological regulations are needed for different areas to accurately assess regional ESVs. The average annual rainfall, transpiration, and evapotranspiration from 2010 to 2014 were used to calculate water yields, with the ratio of mean value of water yield in the study area (Figs. 2(c) and 3(a)) to the mean value for all of China ( Fig. 4(a) and (b)) used as a correction for regional water supply and hydrological regulation of ecosystem services. The water production correction formulae were defined as follows: where WY i denotes the average annual water yield in the study area and WY denotes the average annual water yield nationwide; R i and R a denote the average annual rainfall in the study area and nationwide between 2010 and 2014, respectively; ET i and ET a denote the average transpiration and evapotranspiration for the study area and nationwide from 2010 to 2014, respectively.

NPP spatiotemporal correction factor
Net plant primary productivity (NPP) is an important component of the surface carbon cycle and directly reflects the productive capacity of vegetation communities under natural environmental conditions while also characterizing the quality status of terrestrial ecosystems. In addition, NPP is also a major factor that determines ecosystem carbon sources/sinks and regulates ecological processes, thereby playing important roles in global change and carbon balances (Peng et al. 2016). Hence, NPP differences reflect ESV differences owing to variation among regions and land use types. In this study, the mean NPP value for the study area ( Fig. 3(b)) from 2000 to 2010 was used as a ratio to the mean value for all of China (Fig. 4(c)) to correct for regional gas regulation and nutrient cycling ecosystem services. The NPP correction formula was defined as follows: where N i denotes the average annual NPP for the study area and N denotes the nationwide average annual NPP.

Habitat quality spatiotemporal correction factor
A key goal for sustaining biodiversity is the effective conservation of species' habitats (Rew et al. 2020 provide animals with basic environmental conditions like adequate water, food, shelter, and breeding sites. However, differences exist in the conditions provided for animals in different habitat types. For example, constructed areas provide little habitat for animals and can negatively affect their survival. In addition, the survival and reproductive success of a species partially depends on the number, condition, and spatial arrangement of habitat components (Hirzel and Lay 2008). Further, some animals have migratory habits, wherein the type of land use and distance of the pathway vary during migration, and the resistance values encountered along the way vary in magnitude. A large resistance value refers to poor habitat quality and is not conducive to promoting biodiversity in an area. In this study, ArcGis cost distance was used to calculate the level of resistance of an area to biological migration. Based on previous analyses , forest land was set as the most suitable habitat, while the rest of the land use types were set as cost rasters. The resistance values of land use types (Appendix) were then set as previously  . The habitat quality correction formula was defined as follows: where CD i denotes the degree of migration resistance in this study area (Fig. 3 (c)) and CD denotes the nationwide average degree of migration resistance (Fig. 4(d)).

Reachability spatiotemporal correction factor
Cultural ecosystem services are defined as "the nonmaterial benefits that people derive from ecosystems through spiritual satisfaction, cognitive development, reflection, recreation, and esthetic experiences" (Scholte et al. 2015). Although cultural ecosystems are intangible services, they play important roles in bridging ecosystems and social systems. The Recreation Opportunity Spectrum (ROS) theory posits that the function of cultural services primarily depends on the Recreation Potential Index (RPI) and the Accessibility of Recreation Sites Index (ARSI). In turn, the RPI depends on the attractiveness and accessibility of the area. Further, accessibility reflects the real level of cultural service functions of the ecosystem. The national average road density (Fig. 4(e)) and the average road density (Fig. 3(d)) for the study region were calculated using ArcGis and their ratios were used to reflect the capacity for cultural supply. That is, denser average road values and more convenient transportation accessibility would lead to a stronger capacity for cultural service supply. The reachability correction formula was defined as follows: where RD i denotes the average road density in the study area and RD denotes the national average road density.

Soil conservation spatiotemporal correction factor
Soil erosion carries away large abundances of soil nutrients and destroys soil structures, leading to gradual deterioration of soil ecosystems that then directly affects soil quality, agricultural productivity, and biodiversity. Further, accelerated soil erosion leads to the loss of the soil's ability to hold water, leading to pollution of water sources and damaging water engineering facilities by increasing sedimentation in water (Wynants et al. 2021). Consequently, soil hydraulic erosion intensity was used in this study to correct for soil conservation ecosystem services in the study area. The ratio of the mean value of soil hydraulic erosion intensity for China (Fig. 4(f)) to the mean value of soil hydraulic erosion intensity for a given area (Fig. 3(e)) was used as the regional spatiotemporal correction factor for soil erosion. The soil conservation correction formula was defined as follows: where S i denotes the average erosion modulus for the study area and S denotes the nationwide average erosion modulus.
The basic framework for the study design is shown in Fig. 5.

Land use simulations
To verify simulation accuracy, land use was predicted for 2015 using the land use in 2000 and 2005, revealing Kappa values of 80.2% and 82.6% (Table 1), respectively. The average forecast rates for cultivated lands, forest lands, grasslands, water bodies, built-up areas, and unused lands in 2015 were 78.9%, 79.8%, 83.4%, 85.6%, 75.3%, and 80.7%, respectively. Given that these values were both > 0.75, the simulation results were appropriate and could accurately reflect future land use scenarios. Consequently, the CA-Markov model was used to simulate the 2025 and 2030 (Fig. 6) land use scenarios for the Loess Plateau.

Analysis of land use structures
Grasslands accounted for the largest proportion of land cover in the Loess Plateau (41.7%), followed by cultivated lands (32.3%), forest lands (14.9%), unused lands (6.4%), built-up areas (3.2%), and water bodies (1.4%) (Fig. 7). Among grasslands, moderate grasslands (46.9%) accounted for the largest share (Fig. 8). In addition, dry land primarily constituted the cultivated lands within the Loess Plateau (97.3%). The area of forest land was dominated by forest and shrub lands, with the sum of the two accounting for about 78%. Among built-up areas, rural settlements comprised the largest portion (67.6%), while sandy land comprised the largest portion of unused land (71.6%), and swampland only accounted for 2.4% of unused land. Water body area primarily comprised streams, rivers, and bottomlands (74.9%).

Analysis of land use dynamics
The rate of land use type changes exhibited varying characteristics among different stages during the study period (Table 2), with the most prominent changes being those for built-up areas.
The K value for built-up areas was the highest among land use types (1.11-9.04%) between 1995 and 2025, indicating that built-up area types increased the most over this period. The K value for construction areas was − 0.52% between 2025 and 2030, indicating that built-up areas will begin to decrease during this period. The K value for cultivated land was only positive for the period corresponding to 1995-2000, indicating that cultivated land growth only occurred during this period. Overall, the K values for built-up areas, water bodies, and forest lands were 5.59%, 0.57%, and 0.39% throughout the study period, respectively, indicating increased land uses for these land types. In contrast, the K values for unused land, grassland, and cultivated land were − 0.30%, − 0.25%, and − 0.19%, respectively, indicating decreased land use for these types.

Land use changes and transformation patterns
Built-up areas increased from 1995 to 2018 (Fig. 7), exhibiting a net increase of 131.06 × 10 4 ha (Table 3) and this primarily occurred due to the conversion of cultivated lands and grasslands, with change values of 74.19 × 10 4 ha and 39.84 × 10 4 ha, respectively. Forest lands generally increased (net increase of 1,03.64 × 10 4 ha or 11.7%), and this primarily occurred due to the conversion of grasslands (87.98 × 10 4 ha). Net increases in watershed areas comprised 5.09 × 10 4 ha, primarily from the conversion of grassland (4.29 × 10 4 ha) and unused land (2.62 × 10 4 ha). The net increase in builtup areas was 112.15 × 10 4 ha between 2018 and 2030, representing an increase of 26.9%, with a net conversion of 54.61 × 10 4 ha and 52.21 × 10 4 ha to cultivated land and grassland, respectively. The net increases in forest lands and water bodies were smaller, at 9.17 × 10 4 and 11.8 × 10 4 ha, respectively. The increase in built-up areas, forest lands, and water bodies during these two study periods were primarily concentrated in the core metropolitan area (Fig. 9) within the southeastern region of the Loess Plateau and around the Yellow River, respectively.

Changes in the values of ecosystem services and their contributions
Indicators and data used to correct ESV equivalent weight factor table (Table 4). Land use reclassification data were used for ESV accounting (Table 5).

Changes in ESVs between 1995 and 2030
The high-ESV areas of the Loess Plateau are mostly concentrated in the southeastern region of the research area, around the Yellow River and Qinling Mountains ( Fig. 10(a)). The ESV for the Loess Plateau steadily declines extending from the southeast to the northwest. ESVs exhibited overall increases during 1995-2030 ( Fig. 11(a)), with a total increase of 1,246.83 × 10 8 yuan. In phase 1 (1995-2018), total ESVs increased by 933.97 × 10 8 yuan, with the largest increase of 532.37 × 10 8 yuan occurring between 1995 and 2000. The areas with the most pronounced increases in ESVs were primarily concentrated along the Yellow and Wei rivers ( Fig. 10(b)), in addition to the southeast region of the research area. The areas with the greatest losses in ESVs were primarily concentrated in the northwestern portion of the Loess Plateau and the central region of Xi'an. In 1995, the total ESV was 21,118.49 × 10 8 yuan ( Fig. 11(a)), with regulatory services comprising the largest portion (14,170.50 × 10 8 yuan; 67.1% of the total), followed by supporting services (14.9%), provisioning services (10.6%), and cultural services (7.4%). The total ESV in 2018 was 22,052.46 × 10 8 yuan and this primarily derived from an increase in regulatory services (753.88 × 10 8 yuan) and increases support, cultural, and provisioning services of 84.56 × 10 8 , 54.84 × 10 8 , and 40.69 × 10 8 yuan, respectively (Appendix). In phase 2 (2018-2030), the total ESV increased by 312.86 × 10 8 yuan, while the ESV decreased by 269.75 × 10 8 yuan in 2018-2030. The Hetao irrigation area of Inner Mongolia was estimated to have the most noticeable increase in ESVs from 2018 to 2030, while the Loess Plateau's northern region and core metropolitan area are estimated to see the most significant decreases in  (Fig. 10(c)). In phase 2, support services and cultural services decreased by approximately 46.17 × 10 8 and 9.85 × 10 8 yuan, respectively. In contrast, the values of regulatory and provisioning services increased by 366.18 × 10 8 and 2.71 × 10 8 yuan, respectively. As indicated above, the regions with the most significant changes in ecosystem services were primarily located in the southeastern areas of the Loess Plateau and around rivers.

Contribution of different land use types to ESVs.
Different land use types differentially contributed to ESVs (Fig. 11(b)). The five land use types with the highest contributions to total ESVs throughout the study period were moderate grasslands, water systems, shrublands, scrublands, and forestlands that accounted for 21.0%, 18.1%, 13.8%, 13.0%, and 12.0% of total ESVs, respectively. Although the area covered by water systems was only 0.8% of the study area, the total ESV contributed by water systems was 18.1%. Throughout the study period, the total ESVs increased for water system, forest, shrub, and bare soil land ecosystems, at levels of 699.33 × 10 8 , 941.43 × 10 8 , 467.51 × 10 8 , and 19.55 × 10 8 yuan, respectively. In contrast, the greatest decreases were observed for sparse grass, moderate grassland, and dry land ecosystem services, by 400.6 × 10 8 , 267.84 × 10 8 , and 129.34 × 10 8 yuan, respectively.

Analysis of land use simulation and change
Land use data were evaluated in this study for the Loess Plateau, with land cover datasets deriving from recognized institutions or official statistic databases at a spatial resolution of 30 m and with an accuracy of over 90%. In the land use simulations of this study, a grid size of 300 m × 300 m (converted from a spatial resolution of 30 m × 30 m) was chosen and nine factors that affect land use change were evaluated, followed by the construction of a multi-criteria evaluation method (Appendix) to simulate predicted land use scenarios on the Loess Plateau. The accuracy of the CA-Markov model varied for different land use simulations (Table 1). For example, built-up areas exhibited relatively low prediction accuracy, primarily due to the close association of built-up areas with social and economic development (Cai and Wang 2020). Social and economic development in China is influenced by national policies, and thus, policy changes significantly impact the size and location of construction lands. Kappa tests revealed that longer forecast series led to lower forecast accuracy. For example, the accuracy of the 10-and 15-year forecasts were 82.6% and 80.2%, respectively. Therefore, future forecasts should consider policy influences and the impact of forecast time series on CA-Markov model results as much as possible.
During the rapid economic development and urbanization that has occurred in the past few decades, land use in the Loess Plateau has undergone unprecedented changes (Feng et al. 2016;Fu et al. 2017;Xiu et al. 2021). Croplands, grasslands, and unused lands decreased between 1995 and 2030, while built-up areas, water area, and wetlands increased over this time (Fig. 7). Moreover, forest lands and built-up areas increased by 50.69% and 159.41%, respectively, primarily due to conversion from grasslands and croplands. The reasons underlying these LUCCs are primarily related to policies promulgated by the Chinese government in the region, including returning farmlands to forests and grasslands (Wang et al. 2019), banning free grazing, the Beijing-Tianjin Wind and Sand Source Control Project (Shan et al. 2015;Zhao et al. 2020), and the China Wetland Conservation Initiative (Meng et al. 2017). These ecological conservation measures have altered land use patterns in the Loess Plateau, and the regional ecosystem has gradually developed in a healthy and stable direction. Spatial distribution analysis (Fig. 9) indicated that increases in forest land and built-up Table 4 Indicators and data used to calculate the correction factor  areas were primarily concentrated in the southern region of the Shanxi Province and the core metropolitan area, while increases in arable land, water area, and wetlands were primarily concentrated around the Yellow River. Agricultural intensification is one of the main drivers of land use change and people have reclaimed high-quality farmland close to water sources to generate higher agricultural yields and increase investments in water conservation facilities (Liu et al. 2013). Thus, farmland and water areas have gradually expanded around the river. The accelerated urbanization rate in China after 2000 has resulted in abundant built-up areas encroaching on arable lands. The increases in built-up areas are primarily concentrated in the core metropolitan areas, also indicating that urbanization and infrastructure development in the study area are also accelerating.

Comparison with ESV changes identified in other studies
The ESV valuations of this study are much higher than those based on the amount of value per unit area for two reasons. First, the price basis for expressing the value in this study differs from other studies. Xie et al. (2015) previously used the average net profit per hectare metric for three major crops to develop the ES equivalent factor table. Sharp rises in grain production costs over time, coupled with the Chinese government's moderate increase in macro-control and the associated grain price subsidy policies have led to successive declines or even negative average net profits per hectare for the three major crops. Thus, it is clearly not reasonable to use the average net profit per hectare value for the three main crops across a year to calculate ESV equivalent weight factor tables. Rather, the total value of the three main crops has remained relatively stable and it is consequently more reasonable to use this value to correct the ESV table. Second, differences in the conditions of ecological service functions among different regions are considered, and the strength of different ecosystem service functions is influenced by different ecological processes and conditions (Costanza et al. 1997). Importantly, this study corrected the spatiotemporal heterogeneity of ESVs using six key ecological indicators (Table 4) including the economy, water production, NPP, habitat quality, accessibility, and soil conservation, ultimately obtaining a table of ESV equivalent weighting coefficients that is applicable to the Loess Plateau (Table 6). Furthermore, some issues were identified with the

Table 5
Land use reclassification data for the Loess Plateau The original land use classifications data comes from Liu et al. (2005). Bare soils 65 NPP, precipitation, and evapotranspiration datasets due to them not being the most recent data, temporal inconsistency, poor spatial resolution, and poor accuracy, which arise from their poor accessibility owing to time and effort constraints. Consequently, there is still room for further improvement in these areas. Fortunately, the values of these parameters used herein are independent of each other for correcting the ecosystem service value tables. Thus, their temporal inconsistencies do not affect the corrected results.
ESVs in China have boomed after 2000, but assessment methods primarily include two techniques that are based on equivalent factors or physical quantities, with the former achieving rapid accounting of ESVs due to simple operation and easy comparison of results. Consequently, studies have primarily adopted the equivalent factor method for assessing ESVs. The latest assessment results (Table 6) obtained based on the above studies in China and elsewhere lead to lower values of ecosystem services of the Loess Plateau than are calculated using this table compared to the ESV tables developed by Xie et al. (2015). For example, the value of ecosystem services of the Loess Plateau calculated in this paper for 2018 is 22,052.5 × 10 8 yuan ( Fig. 11(a)), while the value calculated using the table developed by Xie et al. (2015) is 1,516.3 × 10 8 yuan higher than the above value. In particular, soil conservation, gas regulation, and maintenance of nutrient cycling ecosystem service functions are 1,921.87 × 10 8 yuan, 1,246.15 × 10 8 yuan, and 139.76 billion yuan higher in the previous study than in this study, respectively, while esthetic landscape ecosystem service function values are 611.32 × 10 8 yuan lower than in this study. Further examination of these differences revealed that the high soil erosion modulus, low vegetation cover, and high population in the study area compared to the national averages resulted in a smaller accounting value for soil conservation, gas regulation, and the maintenance of nutrient cycling ecosystem service functions (Table 4). The ecosystem value in this study was 1,013.4 × 10 8 yuan higher than that calculated by Jiang et al. (2016). However, the calculation method used by the latter study was subjective and did not consider the spatial and temporal variability of ESVs, thereby leading to low assessment results. Previous studies (Song and Deng 2017;Chen et al. 2020) have assessed ESVs using assessment coefficients primarily based on national averages or global averages, but failed to consider the influence of ecological processes and conditions between different regions. Therefore, this study better reflects the actual ecosystem service values and strength of service functions in the study area by considering regional spatial and temporal heterogeneity. In addition, human perspectives on ecosystem services vary in different socio-economic contexts (Wang et al. 2019). Consequently, this study considers a correction of factors on the "supply side" of ecosystem services that does not reflect human preferences. When using the equivalent factor method in the future, socio-economic factors can be considered and corrections can be made from the "demand side" of ecosystem services. Thus, when using the equivalence factor method in the future, we can consider socioeconomic factors and make corrections from the "demand side" of ecological services.

Impacts of LUCC conversion on ESV changes
Different ecosystems exhibit differences in structure and function (Costanza et al. 1997). For example, when considering the same study area in this investigation, large differences in ESV were generated by dense, moderate, and sparse grasslands of the same region. In this study, secondary land use reclassification (Table 1) was used to more accurately reflect and reliably estimate ESVs. The Loess Plateau's ESVs generally decrease from the northwest to the southeast, mostly following the distribution of forests, thereby demonstrating the significance of woodlands for the ecological security of the region. The increases in ESVs from 1995 to 2030 are consistent with the expansion of forests and watersheds, indicating that these two factors are primarily responsible for observed shifts in the Loess Plateau's ecosystem service values. Decreased ESVs are commensurate with how the expansion of the central urban region was distributed during development. Thus, rapid urbanization has resulted in the land conversion of regions that were once arable or grasslands to built-up areas, thereby leading to severe detrimental effects on the ecology of the region. Notably, the ESVs for the northwestern region of the study area continually decreased during the periods considered in this study, suggesting continual deterioration of that region's ecological conditions. Thus, future research should especially consider the Loess Plateau's ecological capital imbalance. In many developing countries, economic benefits are the primary consideration for land resource management, while ecological benefits are often neglected, thereby increasing pressure on natural resources and ecosystems, leading to threatened regional ecological security (Horlings and Barsden 2011). The ESVs in the Loess Plateau region primarily derive from grasslands, forest lands, and farmlands ( Fig. 11(b)). Across the study period, the value-addition of ecosystem services in the region primarily originated from increases in forest, scrubland, and water system areas that contributed 37%, 18.3%, and 43.9% to the value-added total for ecosystems, respectively. The value of regulatory services increased the most throughout the study period, primarily due to increased area of forest lands and water systems. The contribution of individual ecosystem service functions to total ESVs was highest for hydrological regulation and lowest for maintaining nutrient cycles ( Fig. 12(b)), suggesting that vegetation cover in the area was low and should be allowed to continue increasing to enhance ESVs.
Water systems and wetlands provide the highest ESV per unit area of land use type and play a key role in water supply, hydrological regulation, and cultural landscapes. Thus, their protection can be effective for preventing ESV declines. Food supplies are mainly generated from cultivated lands that play important roles in food production, raw material production, gas regulation, soil formation, soil retention, and nutrient cycling. In addition, forest lands and grasslands play important roles across almost the entire range of ecosystem services. Overall, the conversion of forest lands, grasslands, watersheds, and wetlands to built-up areas and cultivated land results in significant losses of ecosystem services. Consequently, land for construction and cultivation should be rigorously planned to prevent damage to environments and loss of ESVs.

Conclusions
In this study, land use was simulated for the Chinese Loess Plateau in 2025 and 2030 using CA-Markov modeling. A revised framework for ESV assessment was then established to account for modeled ESVs. The results show that (1) built-up and forest land areas increased to the greatest extent between 1995 and 2030. Built-up areas were primarily derived from conversion of cultivated lands and grasslands concentrated in core urban areas, while forest lands increased due to the conversion of grasslands concentrated in the eastern and southern regions of the study area. (2) ESVs generally gradually decreased from the southeastern to the northwestern regions of the study area. Moderate grasslands contributed the largest proportions to ESVs, while the total contributions of water systems, forests, and shrub lands to ESVs gradually increased. (3) Ecosystem services like hydrological regulation, climate regulation, water supply, environmental purification, recreation and culture, and biodiversity increased from 1995 to 2030, while food production ecosystem services significantly declined over that period. (4) ESVs continuously increased throughout the study period, indicating that the ecosystems of the study area are evolving toward healthy and stable states. However, the imbalanced growth and distribution of spatial ecological capital among areas suggests that future development policies must consider the imbalance in ecological capital between regions to ensure sustainable ecological development.