Evaluating the influences of urban expansion on multiple ecosystem services in drylands

Effectively estimating the influences of urban expansion on multiple ecosystem services (ESs) is of great importance for improving urban planning in drylands. However, there are some shortcomings in the existing urban expansion models, which lead to great uncertainties in the assessment of the influences of urban expansion on multiple ESs. This study sought to effectively estimate the influences of urban expansion on multiple ESs in drylands. We combined the land use scenario dynamics-urban model and ES models to estimate the influences of urban expansion on habitat quality, food production, water retention, air quality regulation, sand fixation, and aesthetics. The results showed that (1) our method can effectively evaluate the influences of urban expansion on multiple ESs in drylands, and the accuracy increased by more than 20% on average. (2) Under the impacts of future urban expansion, the annual average loss per unit area of habitat quality, food production, water retention, air quality regulation, sand fixation, and aesthetics will increase. (3) Future urban expansion will occupy more cropland and grassland, which will be the main reason for the intensified degradation of multiple ESs. We suggest that the local government should focus on protecting key ecosystems, especially grasslands, from urban expansion to mitigate the potential influences on multiple ESs in drylands.


Introduction
Ecosystem services (ESs) are the benefits that humans obtain from ecosystems (MEA, 2005). The loss of multiple ESs caused by urban expansion refers to the decrease of multiple ESs provided by a certain ecosystem when it is occupied by urban land and its structure and functions are altered during urban expansion (Gong and Liu 2021;Wang et al. 2021). For example, after urban expansion occupies grassland, it will lead to the loss of water retention service, air quality regulation service and sand fixation service provided by grassland. Previous studies have revealed that urban expansion has led to the loss of regulating services, provisioning services, and supporting services (Gomes et al. 2020;Zhang et al. 2021;Hou et al. 2021), which has become an obstacle to regional sustainable development (MEA 2005;Nelson et al. 2010). Drylands refer to areas with water shortages as the main feature, in which productivity and nutrient circulation are restricted by the water supply. Drylands account for 41.3% of the world's land area, in which 38% of the world's population lives (MEA 2005;Cherlet et al. 2018). Drylands ecosystems have not only largely contributed to the global land carbon sink but also contain approximately one-third of biodiversity hotspots and provide habitat for 28% of endangered species (Fu et al. 2021;Stafford-Smith and Metternicht 2021). With the rapid expansion of urban land, the reduction in ESs is aggravated in drylands (Reynolds et al. 2007;Xie et al. 2018;Adelisardou et al. 2021). In the future, drylands will still undergo large-scale urbanization (Seto et al. 2012;Zhou et al. 2019). Therefore, the effective assessment of the potential influences of future urban expansion on multiple ESs in drylands is of great importance for improving urban planning and promoting regional sustainable development.
At present, researchers have analyzed the influences of urban expansion on multiple ESs in drylands. For example, Adelisardou et al. (2021) found that carbon storage and water yield showed a downward trend as the urban land area increased in Jiroft plain, Iran from 1996 to 2016. Narducci et al. (2019) investigated the impacts of urbanization on multiple ESs in the Boise metropolitan area, USA, and found that urbanization has significant negative impacts on food production, water quality, air quality, and habitat for species. Lyu et al. (2018) evaluated the effects of urbanization on ESs during the period of 1989-2015 in the northern part of Ningxia, China, and the results showed that urban expansion exerted negative impacts on carbon storage, crop production, sand fixation, habitat quality, and nutrient retention. Xie et al. (2018) assessed the impacts of urban expansion on water conservation, food production, air quality regulation, habitat quality and carbon storage in Beijing, China, during 2013-2040 and found that urban expansion will lead to the simultaneous loss of multiple ESs in the future.
However, there are still some uncertainties in existing studies. This is mainly because the existing urban expansion models have certain deficiencies in simulating future urban expansion. For example, the urban expansion model used by Xie et al. (2018) is based on the traditional cellular automation (CA) model. First, the traditional CA model assumes that the way of urban expansion remains unchanged from the past to the future and uses static model parameters in the simulation. Second, errors from data sources (such as field survey errors, digitization errors, and data conversion errors) will transmit and accumulate constantly in the CA model, thereby increasing the uncertainty associated with simulating future urban expansion (Zhang et al. 2015).
Improving the existing urban expansion model with the data assimilation method can lead to the more accurate simulation of future urban expansion, thereby reducing the uncertainties of the evaluation results of future urban expansion on the loss of multiple ESs. The data assimilation method directly or indirectly integrates multiple sources and multiple resolution observations by means of observation operators while weighing the uncertainties in the model and the observations to minimize the error of the entire system (Sakov and Bertino 2011;Li et al. 2020). This method has been widely used in earth system sciences (Li et al. 2020). The ensemble Kalman filter (EnKF) is a widely used data assimilation method (Li et al. 2020). At present, previous studies have tried to use the EnKF method to improve the urban expansion model. For example, Zhang et al. (2015) optimized the CA urban expansion model using the EnKF method and found that the EnKF method can reduce error transmission and dynamically optimize the relevant parameters of the CA model, thereby improving the simulation capability of the urban expansion model. However, there is still a lack of studies using an improved urban expansion model based on data assimilation to evaluate the influences of future urban expansion on multiple ESs.
The objective of this paper is to improve the land use scenario dynamics-urban (LUSD-urban) model with data assimilation and then combine it with ES models to evaluate the potential influences of future urban expansion on multiple ESs. First, the Hohhot-Baotou-Ordos-Yulin (HBOY) urban agglomeration, China was chosen as the study area and we mapped six key ESs by employing ES models. Then, we used the EnKF method to improve the LUSD-urban model and simulated future urban expansion. Finally, we estimated the potential influences of future urban expansion on multiple ESs. This study will provide an effective method for assessing the impacts of urban landscape dynamics on the environment, which is of great significance for monitoring urban landscape and evaluating their sustainability (Wu 2014;Liu and Wu 2022).

Study area
The HBOY urban agglomeration is located in the middle of the drylands of northern China (E 106.5-122.3°, N 36.8-42.7°), with a total area of 1.75 × 10 5 km 2 (Fig. 1). The average elevation in the study area is approximately 1300 m. The climate type is a temperate continental monsoon climate, with a 45-year average temperature of approximately 7.48 ℃ and a 45-year average precipitation of approximately 349.46 mm (Song et al. 2020). There are 30 cities and 439 towns in the HBOY urban agglomeration.
In recent decades, the HBOY urban agglomeration has experienced rapid urban expansion, which has caused serious impacts to the regional ESs, thereby bringing new challenges to biodiversity conservation and sustainable development (He et al. 2014;Sun et al. 2017;Song et al. 2020). For example, He et al. (2014) showed that urban expansion in the HBOY region from 1992 to 2014 caused the loss of 1178 km 2 of natural habitat, posing negative impacts to regional species. Sun et al. (2017) found that regional urban expansion from 1990 to 2013 caused the four ESs of habitat quality, food production, meat production, and carbon sequestration to decline. Song et al. (2020) showed that urban land area increased from 151.29 km 2 in 1990 to 1230.86 km 2 in 2017 in the HBOY area, which caused the natural habitat quality to decline. With the development of the regional economy and the relevant national policies (such as the plan for HBOY urban agglomeration and the "One Belt and One Road" initiative), the HBOY will also undergo rapid urban expansion in future, which will pose a huge threat to biodiversity and the surrounding ecosystem (Song et al. 2020;National Development and Reform Commission, 2018). Therefore, the HBOY urban agglomeration provides an appropriate study area for evaluating the potential influences of future urban expansion on multiple ESs in drylands.

Data
In this study, urban land data, land use/land cover (LULC) data, precipitation data, socioeconomic data, and basic geographic information data were used. Among them, the LULC data for 1990 comes from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (www. resdc. cn/ data. aspx? DATAID= 283). The resolution of LULC data is 30 m, and the accuracy is above 90% (Liu et al. 2010). The urban land data during the period of 1990 to 2017 come from the dataset for the HBOY region obtained by the Google Earth Engine. The dataset includes urban land in 1990, 2000, 2010, and 2017, respectively. The resolution of urban land data is 30 m, and the average Kappa is above 0.80 (Song et al. 2020). The annual precipitation data come from the spatial interpolation dataset of annual precipitation in China since 1980 from the Resource and Environmental Science Data Center of the Chinese Academy of Sciences (http:// www. resdc. cn/ data. aspx? DATAID= 229). Socioeconomic data (such as GDP and total population) come from the Statistical Yearbook of Hohhot, Baotou, Ordos, and Yulin in 2017. The administrative boundaries, highways, railways, and national roads in the HBOY region come from the National Basic Geographic Information Center (http:// ngcc. sbsm. gov. cn/). To ensure data consistency, all data used a unified Albers projection and were resampled to 100 m.

Simulating urban expansion from 2017 to 2050
Improved LUSD-urban model based on EnKF LUSD-urban model The LUSD-urban model consists of two parts: the urban land demand module and the urban land space allocation module. In the demand module, a linear regression model is established based on the urban land area and urban population to predict the total urban land demand in the future. In the space allocation module, whether nonurban grids are converted into urban grids is mainly affected by their suitability, the inheritance attributes of the LULC class, neighborhood effects, and random perturbations in urban expansion. Details of the model formulas can be found in He et al. (2006He et al. ( , 2016. EnKF model The EnKF model can be divided into two stages: prediction and update. In the prediction stage, the simulation of the process model is performed using the initial state set obtained by random sampling to obtain the prediction set, and the Kalman gain is calculated by employing the prediction set. In the update stage, the prediction set is updated by the observations and Kalman gain, and the updated result is used as the analysis value set. The mean of the analysis value set is the posterior estimate of the model and will continue to act as the new state set that participates in the next cycle (Sakov and Bertino 2011;Zhang et al. 2015). Please refer to Online Appendix A for specific calculation formulas.
Combining the LUSD-urban model and the EnKF model The basic principle of the coupling model is to input the simulation results of the urban expansion model and observation data into the EnKF model for data reprocessing and then use observation data to continuously optimize the simulation results and related parameters through the Bayesian method, thereby realizing the dynamic update of model parameters and improving model accuracy. The main advantage of the coupled model is that it only needs to reanalyze existing data to realize improvement of the model simulation accuracy.
First, we used the factor data (i.e., distances to city centers, distances to highways, distances to rivers, and distances to national roads), historical urban expansion data and the LUSD-urban model to build the urban expansion model and obtain the initial parameters (Fig. 2).
Second, we put the initial parameters into urban expansion model (LUSD-urban model) to get the urban land simulation results. The EnKF model requires the input data to be non-Boolean, while the simulation results of the LUSD-urban model only have two states, namely, urban (pixel state is 1) and non-urban (pixel state is 0). According to Zhang et al. (2015), we divided the study area into several grids, and each grid included a certain number of pixels. We calculated the urban expansion intensity of each grid as follows: (1) where t i,j is the urban expansion intensity for grid (i,j) at time t; ∑ m×m con � s i,j = urban � is the total number of urban pixels in grid (i,j); and m × m is the total number of pixels in grid (i,j).
Third, the simulation results and parameters (the weights of factors) of LUSD-urban were corrected by integrating the EnKF model and the actual urban land data. The main steps were as follows: (a) Several grids are randomly selected as observation points, and the urban expansion intensity of those grids is the observation value. (b) The LUSD-urban model is used to simulate urban land in the same year of data assimilation, and the urban expansion intensity of the simulation results is calculated as the prediction value. (c) The Kalman gain is calculated using the observation value and the prediction value. (d) The simulation results and the parameters of the LUSD-urban model are corrected using the Kalman gain and observations. The specific formula can be found in the Online Appendix B.
Finally, we put the corrected results and parameters into the LUSD-urban model as the initial inputs of the next cycle until all observations have been assimilated. After assimilation, the improved LUSD-urban model can be used to simulate future urban expansion.

Simulation accuracy evaluation
In order to verify the validity of the model, we used Kappa coefficient, location Kappa coefficient, and figure of merit (FoM) to verify the model simulation results. Among them, the Kappa coefficient evaluates the accuracy of the simulation results by quantifying the overall correlation between predicted data and actual data (van Vliet et al. 2011). The formula is as follows: where P o is the probability that the simulated data agree with the actual data; P e is the probability that the simulated data agree with the actual data by accident. The Kappa coefficient has a range of − 1 to + 1. The larger value means that the higher the consistency between the simulation result and real data, and the higher the simulation accuracy. Generally, Kappa > 0.8 indicates very high accuracy, 0.6 < Kappa < 0.8 indicates good accuracy, 0.4 < Kappa < 0.6 indicates poor accuracy, Kappa < 0.4 indicates very poor accuracy ). The formula for P o and P e are as follows, where A is the number of pixels that are actually urban pixels and simulated as urban pixels; B is the number of pixels that are actually urban pixels but simulated as non-urban pixels; C is the number of pixels that are actually non-urban pixels but simulated as urban pixels; D is the number of pixels that are actually non-urban pixels and simulated as non-urban pixels. The location Kappa coefficient has a range of − 1 to + 1, and a larger value indicates the consistency of spatial location between the simulation result and real data is higher.
According to Chen et al. (2016), the formula for FoM is as follows, where Hits is the number of pixels that are observed change and simulated as change; Misses is the number of pixels that are observed change but simulated as non-change; FalseAlarms is the number of pixels that are observed non-change but simulated as change. The FoM has a range of − 1 to + 1, and a larger value indicates a higher simulation accuracy.
In our study, we put LULC in 2000, elevation, slope, distance to the city centers, distance to rivers, distance to highways, distance to railways, and distance to national roads into the LUSD-urban model and used the Monte Carlo method to obtain the initial parameters, including the weights of factors and the simulation result of urban land in 2010. Then, the simulated urban land in 2010 and actual urban land in 2010 were transformed into urban expansion intensity and put into the EnKF model for data assimilation. Finally, the assimilated simulation results and parameters were input into the LUSD-urban model to predict the urban land area in 2017, and the validation of the model was verified with the actual urban land area in 2017. The results showed that Kappa, Kappa location and FoM were 0.64, 0.64 and 0.19, respectively, indicating that the model can effectively simulate future urban land expansion (Online Appendix B).

Projecting urban expansion from 2017 to 2050
According to He et al. (2016), the regression model was used to forecast urban land demand in the future. In this study, we used the urban land area and urban population in 1990, 2000, 2010, and 2017 to build a linear regression model. The formula is as follows: where S is the urban land area; Pop is the urban population. The R 2 value of the linear regression model is 0.97.
Then, we used urban population and time to construct a regression model to forecast the urban population in 2050. The formula is as follows: where Pop is the urban population and y is the year. The R 2 value is 0.98.
Finally, we predicted that the urban land area in the HBOY region will reach 2498.25 km 2 in 2050 based on Eqs. (16) and (17). On this basis, the improved LUSD-urban model was used to project the urban expansion in the HBOY urban agglomeration from 2017 to 2050.

Mapping ESs
According to the availability of data and the Millennium Ecosystem Assessment, we chose six key ESs, including supporting service (habitat quality), provisioning service (food production), regulating service (water retention, air quality regulation and sand fixation), and cultural service (aesthetics).

Habitat quality
The Habitat Quality Module of the InVEST (the Integrated Valuation of Ecosystem Services and Tradeoff) model was used to quantify the habitat quality. This module calculates habitat quality based on the external influence factors and the sensitivity of a given LULC type. The main formula (Sharp et al. 2015) is as follows: where HQ ij is the habitat quality for pixel j in habitat type i; H i is the suitability in habitat type i; D z ij is the total threat for pixel j in habitat type i; K is the half-saturation constant and Z is the scaling parameter. The related parameters are obtained from Song et al. (2020).

Food production
Referring to Li et al. (2014), we mapped the food production by LULC data. In the method, cropland yields vegetables, oil and grain, and grassland yields milk and mutton. The formula can be expressed as follows: where FP ij is the total food production service for grid j under LULC type i; A j is the area of grid j; and p fi is the per unit area yield for food f in LULC type i.
The p fi parameters can be found in the Supplementary file (Online Appendix C).

Water retention
Referring to Yang et al. (2015), the capability of vegetation (such as cropland, grassland, and forest land) to intercept surface runoff was calculated to quantify water retention service. The formula is as follows: where WR ij is the total water retention for grid j in LULC type i; A j is the area of grid j; P j is the 27-year (1990-2017) averaged precipitation for grid j. K is the surface runoff coefficient, which is set to 0.6 based on Yang et al. (2015). R ij is the proportion of surface runoff relative to total rainfall. The R kxy parameters can be found in the Supplementary file (Online Appendix D).

Air quality regulation
Referring to Landuyt et al. (2016), the capability of vegetation to adsorb the PM 10 (particulate matter with a particle size below 10 μm) was used to quantify the air quality regulation service. The formula is as follows: where AQR ij is the total retention volume of PM 10 for grid j in LULC type i; A xy is the area of grid j; and PM ij is the retention volume of per unit area PM 10 for grid j in LULC type i. The PM ij parameters can be found in the Supplementary file (Online Appendix D).

Sand fixation
According to Wei et al. (2018) and Zhao et al. (2017), we used the difference between the potential amount and the actual amount of sand erosion by wind to map the sand fixation service. The main formula is as follows: where SF j is the amount of sand fixation service for pixel j; sf jp is the potential amount of sand erosion by wind for pixel j; sf ja is the actual amount of sand erosion by wind for pixel j.

Aesthetics
According to Cui et al. (2019), the visual quality index is used to quantify the aesthetics service, which represents the aesthetic appeal of natural landscapes to tourists. The formula is as follows: where VQI j is the total visual quality index score for pixel j, with the range of 0-5; T j is the terrain factor score for pixel j; W j is the water bodies factor score for pixel j; V j is the natural vegetation factor score for pixel j;U j is the urbanization factor score for pixel j; and L j is the natural landscape accessibility factor score for pixel j.

The loss of multiple ESs during urban expansion
At present, there are two types of methods to quantify the loss of multiple ESs caused by urban expansion, namely direct measurement and indirect calculation. In terms of direct methods, the land use data and processbased models were combined for measurement (Chen et al. 2019). Among indirect methods, the correlation analysis was often used to analyze the relationships between ESs (Xie et al. 2018;Lyu et al. 2021). Considering the availability of data and the simplicity of the (14) SF j = sf jp − sf ja (15) VQI j = T j + W j + V j + U j + L j method, we refer to Xie et al. (2018) and adopt the correlation analysis method to analyze the relationships of ESs during urban expansion. First, the spatial patterns of the six ESs were mapped by the ES models outlined in Sect. 3.1, and their changes were calculated from 1990 to 2017 and 2017-2050. Changes were calculated as follows: where ΔES i,j is the change in ecosystem service i in township j from year t 1 to year t 2 ; ES t 1 i,j is ecosystem service i in township j in year t 1 ; and ES t 1 i,j is ecosystem service i in township j in year t 2 . If ΔES i,j is more than zero, it indicates that ESs are enhanced. If ΔES i,j is negative, it indicates that ESs are degraded and that urban expansion has a negative impact on ESs.
Then, Pearson's correlation coefficient was calculated to analyze the correlation between ESs. The calculation of the correlation coefficient is as follows: where R is the correlation coefficient between ecosystem service X and Y; cov(X, Y) is the covariance between ecosystem service X and Y; D(X) is the variance of ecosystem service X; and D(Y) is the variance of ecosystem service Y.

Urban expansion from 1990 to 2050 in the HBOY Region
The HBOY region experienced rapid urban expansion from 1990 to 2017. The urban land area increased from 151.29 to 1230.86 km 2 , representing a 7.14fold increase, with an average annual growth rate of 19.29%. Among them, Baotou, Hohhot, Ordos, and Yulin are the four cities with the largest increase in urban land, with areas of 320.91 km 2 , 214.77 km 2 , 113.78 km 2 , and 80.50 km 2 , respectively. These cities accounted for 67.62% of the total area of urban expansion in the HBOY region (Fig. 3A).
The HBOY will continue to undergo rapid urban expansion from 2017 to 2050. The area of urban land will reach 2498.25 km 2 in 2050, with an increase of (16) � 102.97%. The new urban land will be mainly distributed in Hohhot, Baotou, Yulin, and Ordos, with areas of 370.38 km 2 , 320.53 km 2 , 114.38 km 2 , and 114.07 km 2 , respectively. These cities will account for 72.54% of the total area of newly increased urban land throughout the whole region (Fig. 3B).
The loss of multiple ESs during urban expansion from 1990 to 2017 Urban expansion from 1990 to 2017 led to a decline in the six ESs in the HBOY. Among them, the loss of food production was the largest, which was reduced from 132.87 to 131.77 t/km 2 , for a total decrease of 0.83%. The loss of aesthetics was the smallest, which was reduced from 3.85 to 3.84, decreasing by 0.26%. Water retention, air quality regulation, sand fixation and habitat quality decreased by 0.54%, 0.43%, 0.67% and 0.39%, respectively (Fig. 4). Spatially, urban expansion in Baotou and Hohhot had the most serious impacts on the evaluated ESs (Fig. 4). Urban expansion in Baotou from 1990 to 2017 led to reductions in habitat quality, food production, water retention, air quality regulation, sand fixation and aesthetics by 8.09%, 14.98%, 9.50%, 7.63%, 7.04% and 5.17%, respectively (Table 1). Hohhot decreased by 4.60%, 12.32%, 7.02%, 4.51%, 5.42% and 4.15%, respectively (Table 1). The three pairs of ESs shown strong positive relationship during urban expansion from 1990 to 2017. The related three pairs of ESs were habitat quality and air quality regulation, food production and water retention, and water retention and aesthetics. Specifically, the correlation coefficients of the three pairs of ESs in 1990 and 2017 were greater than 0.5 and passed the significance level at 0.01 (Fig. 5a,  b). Among all services, the Pearson's correlation coefficient between water retention and aesthetics increased the most with the urban expansion from 2019 to 2017, which increased from 0.57 to 0.62.

The loss of multiple ESs during urban expansion from 2017 to 2050
With the expansion of urban land in the future, the six key ESs will continue to degrade. From 2017 to 2050, food production, water retention, habitat quality, air quality regulation, sand fixation, and aesthetics will be reduced by 1.34%, 0.81%, 0.65%, 0.60%, 0.95%, and 0.35%, respectively. Among them, the loss of food production will be the most serious and decreased from 131.77 to 130.01 t/km 2 . The loss of aesthetics was the smallest and decreased from 3.84 to 3.83 (Fig. 6). Spatially, the influences of urban expansion on regional ESs in Baotou and Hohhot will still be the most serious (Fig. 6). Among them, urban expansion in Hohhot from 2017 to 2050 will lead to reductions in food production, water retention, habitat quality, air quality regulation, sand fixation, and aesthetics by 8.91%, 26.75%, 13.55%, 7.99%, 11.02%, and 6.56%, respectively ( Table 2). Those of Baotou will decrease by 11.52%, 20.65%, 12.04%, 8.94%, 8.73%, and 5.19%, respectively ( Table 2).
The five pairs of ESs shown strong positive relationship during urban expansion from 2017 to 2050. Among them, habitat quality and air quality regulation, food production and water retention, and water retention and aesthetics will still show strong positive relationship. At the same time, food production and aesthetics and water retention and air quality regulation will also reveal strong positive relationship in the future. Specifically, the related ESs pairs will reveal a significant strong positive correlation in 2017 and 2050 (Pearson's correlation coefficients will be greater than 0.5 and pass the significance at 0.01) (Fig. 5b, c). Among all services, the Pearson's correlation coefficient between water retention and air quality regulation will increase the most during the urban expansion from 2017 to 2050, which will increase from 0.54 to 0.59. Urban expansion will aggravate the loss of multiple ESs in the future. With the increase in urban areas, the annual average loss of ESs per unit area in 2017-2050 will be higher than that in 1990-2017. In the whole region, the annual average loss per unit area between 2017 and 2050 of food production, water retention, habitat quality, air quality regulation, sand fixation, and aesthetics will be 4.09 t/km 2 , 2.48 t/km 2 , 0.0198, 1.81 kg/km 2 , 2.88 t/km 2 , and 0.0105, respectively. The annual average loss of ESs per unit area were 3.08 t/km 2 , 2.02 t/km 2 , 0.0145, 1.61 kg/ km 2 , 2.49 t/km 2 , and 0.0099 in 1990-2017, respectively.

Discussion
Combining the improved LUSD-urban model and ES models can effectively assess the potential influences of future urban expansion on multiple ESs in drylands Improving the LUSD-urban model by the EnKF method can simulate urban expansion more effectively, thereby reducing the uncertainties of evaluating the influences of urban expansion on multiple ESs. This is mainly because the EnKF method incorporates more information (e.g., observation data) and constantly optimizes the LUSD-urban model parameters and the simulation results, which reduces the uncertainties introduced by the static model parameters in the LUSD-urban model and eliminates error accumulation and transmission in the model (Li et al. 2020). Moreover, the EnKF method enhances the accuracy of simulation and prediction by reanalyzing observational data and then mining more data information (Hunt et al. 2007;Carrassi et al. 2018), which enables decision-makers to obtain more precise data without requiring extra expertise. Therefore, the improved LUSD-urban model is theoretically better than the original model in simulating future urban expansion.
In order to verify the effectiveness of our method, we compared the accuracy of the improved LUSD-urban model and the original LUSD-urban model in terms of the simulation of urban expansion and the evaluation of the impacts of urban expansion on multiple ESs. According to the simulation of urban expansion, the accuracy of the improved urban expansion model is better than the original LUSD-urban model. In terms of impact evaluation, the evaluation results using the improved LUSD-urban model and the ES models were closer to the true value. Compared with the original LUSD-urban model, the results of combining the improved LUSD-urban and the ES models showed that urban expansion reduced habitat quality, food production, water retention, air quality regulation, sand fixation, and aesthetics by 11.99%, 50.89%, 26.77%, 17.94%, 41.04%, and 11.72%, respectively, which were closer to the influences of urban expansion on ESs from 2010 to 2017 in the HBOY region (Table 3). In addition, according to He et al. (2016), the relative error was also selected to further verify the effectiveness of our method. The outcomes revealed that the relative errors of our method were − 24.64%, 44.94%, 18.66%, 1.47%, 16.06%, and 8.02%, respectively, which were lower than those of the original LUSD-urban model, and the simulation accuracy was improved by more than 20% on average (Table 3). Furthermore, we compared the simulation results of the improved model and the original model with actual LUCC data from 2010 to 2017 and found that the occupancy of other land types by urban land in the simulation results of the improved model was closer to actual data (Table 4). This was the main reason why the improved model was used to evaluate the impacts of urban expansion on multiple ESs with higher accuracy (Fig. 7).
The occupation of cropland and grassland by urban expansion is the main reason for the loss of multiple ESs Urban expansion occupies a large number of cropland and grassland ecosystems, which is the main reason for the loss of multiple ESs in the HBOY urban agglomeration. In the whole region, the area of cropland and grassland occupied by the urban expansion from 1990 to 2017 was 365.22 km 2 and 321.13 km 2 , respectively, accounting for 33.82% and 29.74% (Table 5). We found that the loss of habitat quality, food production, water retention, air quality regulation, sand fixation, and aesthetics services caused by the occupation of cropland and grassland by urban expansion in the HBOY from 1990 to 2017 accounted for more than 70% of the total loss (Table 5). From 2017 to 2050, the area of cropland and grassland occupied by the urban expansion will be 608.53 km 2 and 402.98 km 2 , accounting for 48.02% and 31.80% (Table 5). This will cause ESs to drop by more than 75% (Table 5). In addition, urban expansion in the future will occupy more cropland and grassland, which will be the dominating reason for the loss of the six ESs. We found that compared with the period from 1990 to 2017, the amount of cropland and grassland occupied by urban expansion will increase by 47.38% from 2017 to 2050, which will lead to an average 1.44-fold increase in ES losses. Grassland ecosystems should be especially protected in the process of urbanization in drylands. Our research found that the occupation of grassland ecosystems by urban expansion has more serious impacts on multiple ESs. Compared with other ecosystems, the loss of ESs per unit area caused by the reduction of grassland ecosystems in the HBOY was highest, such as habitat quality, water retention, air quality regulation, and sand fixation services (Table 5).
Other studies have found that urban expansion in drylands mainly occurred in grasslands and caused the loss of multiple ESs, which is consistent with the conclusions of our study. For example, Sun et al. (2017) showed that urban expansion of the Hohhot-Baotou-Ordos area in the drylands of northern China took up 132.32 km 2 of grassland during 1990-2013, accounting for 35.02%, which led to regional meat production and carbon storage dropped by 100% and 46.06%. Song et al. (2020) also showed that the HBOY urban agglomeration will have more than 60% of urban land converted from grassland from 2017 to 2050, which will cause natural habitat quality to drop by more than 75%. Lyu et al. (2018) found that 42.34% of the newly added urban land in the northwestern region of Ningxia Hui Autonomous Region, China from 1989 to 2015 came from grassland, while carbon storage and nutrient retention services had declined. For non-arid areas, urban expansion often    has more serious impacts on the forest ecosystems (Pickard et al. 2016;Xie et al. 2018).

Policy implications
Our study shows that the loss of multiple ESs will be unavoidable in the future, especially for food production and water retention. Urban expansion led to a 0.83% decrease in food production from 1990 to 2017, with the highest impact reaching 14.98% in Baotou, followed by a 0.54% decrease in water retention, with the highest impact reaching 7.02% in Baotou. Urban expansion will lead to a 1.34% reduction in food production during 2017-2050, with the highest impact of 26.75% in Hohhot, followed by a 0.81% decrease in water retention, with the highest impact of 13.55% in Hohhot. In addition, both food production and water retention showed strong positive relationship during these two periods. The large loss of food production will certainly bring great challenges to regional food security, and the large decline in water retention caused by urban expansion will certainly become the main bottleneck of urban development in drylands (Liu et al. 2019). Therefore, for the drylands, the impacts of urbanization on water retention service should be fully considered in future urban planning and management, and urban construction should be strictly based on the regional water resources carrying capacity.
In the future, the occupation of the key ecosystem of grassland by urban land should be strictly controlled in the HBOY, so as to alleviate the impacts of urbanization on water retention service, especially in Hohhot and Baotou. For example, Baotou should avoid the urban expansion that occupies grassland in the north of the city to protect the ecological environment of the Daqingshan area, and the ecological environment along the Yellow River should be maintained during the urban expansion to the south. Hohhot should also avoid new urban land in the east and north of the city to reduce the impacts on the ecological environment of the Daqingshan area, thereby enhancing the water retention capacity of the regional ecosystem. In addition, future efforts should also aim to control the scale and optimize the spatial pattern of cities to cut down the occupation of cropland and grassland ecosystems by urban expansion.

Future perspectives
Our study proposed a new approach to effectively assess the potential impacts of future urban expansion on multiple ESs in drylands. We combined the improved LUSD-urban model with ES models and used the HBOY urban agglomeration to verify the method effectiveness. The results revealed that our method improved the accuracy of the assessment by more than 20% on average. Our model improved the accuracy of urban landscape simulation, which can provide a reliable and powerful tool for urban landscape monitoring and evaluation. Our findings that urban expansion in drylands will exacerbate the loss and trade-offs of multiple ESs, particularly in food production and water retention, enlighten us to pay attention to urban landscape design in rapidly urbanizing drylands to maintain regional sustainable development (Quintas-Soriano et al. 2016;Asadolahi et al. 2018;Adelisardou et al. 2021).
There are also some uncertainties in this study. First a simple linear regression model was used to determine the future urban population and urban land demand. However, the relationships between economy, population, environment, policy and urban expansion are very complex. Second, the quantification of ESs supply based on LULC data in this study was relatively simple, which may lead to the underestimation of the impacts of urban expansion on the studied ESs. For example, the influences of urban land on ESs due to the occupation of water bodies were ignored in this study. However, this does not impact the dependability of the results of our study. This is because these quantitative methods have been fully applied and verified in existing studies (Yang et al. 2015;Landuyt et al. 2016;Xie et al. 2018), and the occupation of water bodies by regional urban expansion is considered to be relatively minimal, so the influences of urban expansion on multiple ESs can be accurately depicted in our study. Third, our model only evaluated the loss of ESs caused by the occupation of other types of land by urban land and did not consider the conversion between other types of land, which will lead to some uncertainties in the evaluation results. For example, the reduction of cultivated land caused by urban expansion can be supplemented by forest land/grassland. In addition, the reasons why urban expansion changes land use and then results in changes in the relationship between ecosystem services are complex, as well as the relationship between ecosystem services is not necessarily linear. The linear correlation analysis used in our paper can only reflect the changes in the correlation between ESs.
In the future, we will use a more mechanistic model (such as system dynamics model) to predict the future urban land demand based on the analysis of the driving mechanism of urban expansion. In addition, we will further explore more effective and accurate methods of mapping ESs, thereby reducing the uncertainties of the assessment results. We will further improve the evaluation model and set up conversion rules between different land types to comprehensively evaluate the impacts of LUCC on multiple ESs. We can also link the urban expansion model with the ecological process-based model (such as integrated biosphere simulator ecosystem model, Biome-BGC) to quantitatively measure the process and speed of the loss of multiple ESs caused by urban expansion.

Conclusions
Our method can be used to effectively assess the potential influences of future urban expansion on multiple ESs in drylands. First, the improved LUSD-urban model was obviously superior to the original LUSD-urban model in the urban land simulation. The Kappa value of the improved LUSDurban model in simulating urban land in 2017 was 0.64, which was higher than the 0.59 of the original model. Second, compared with the original model, the influences of regional urban expansion on multiple ESs evaluated by combining the improved LUSD-urban model and ES models were closer to the actual value, with smaller relative errors. The accuracy of the assessment was improved by more than 20% on average. Among them, the accuracy associated with evaluating the influences of urban expansion on water retention and sand fixation services improved the most at 61.19% and 29.88%, respectively.
Under the influences of urban expansion from 1990 to 2017, food production, water retention, air quality regulation, habitat quality, sand fixation and aesthetics in HBOY showed a downward trend, which decreased by 0.83%, 0.54%, 0.43%, 0.39%, 0.67% and 0.26%, respectively. During 2017-2050, urban expansion will accelerate the decline of these six key ESs, and the annual average loss per unit area will reach 4.09 t/km 2 , 2.48 t/km 2 , 1.81 kg/km 2 , 0.0198, 2.88 t/km 2 , and 0.0105, which will be higher than that of the 3.08 t/km 2 , 2.02 t/km 2 , 1.61 kg/km 2 , 0.0145, 2.49 t/km 2 , and 0.0099 losses that occurred during 1990-2017. Urban expansion in the future will occupy more cropland and grassland, which will be the dominating reason for the loss of multiple ESs. Compared with 1990-2017, the cropland and grassland occupied by urban expansion will increase by 47.38% in 2017-2050, which will lead to an average increase of 1.44 times the loss of ESs.
Therefore, we suggest that the future urban expansion encroaching on large amounts of grassland should be strictly controlled in urban planning and management to reduce the risk of the impacts of urban expansion on multiple regional ESs, especially for cities undergoing rapid urban expansion and serious impacts on ESs, such as Hohhot and Baotou.