Evaluating the In�uences of Urban Expansion on the Concurrent Loss of Multiple Ecosystem Services in Drylands

: 24 Context Effectively estimating the influences of urban expansion on multiple ecosystem services 25 (ESs) is of great importance for improving urban planning in drylands. However, there are some 26 shortcomings in the existing urban expansion models, which lead to great uncertainties in the 27 assessment of the influences of urban expansion on the concurrent loss of multiple ESs. 28 Objectives This study sought to effectively estimate the influences of urban expansion on the 29 concurrent loss of multiple ESs in drylands. 30 Methods We combined the improved the urban expansion model and ES models to estimate the 31 influences of urban expansion on five key ESs, including food production (FP), water retention 32 (WR), air quality regulation (AQR), natural habitat quality (NHQ), and landscape aesthetic (LA). 33 Results The results showed that (1) our method can effectively evaluate the influences of urban 34 expansion on the concurrent loss of multiple ESs in drylands, and the accuracy increased by more 35 than 20% on average. (2) Under the effect of future urban expansion, FP, WR, AQR, NHQ and LA 36 will accelerate the decline. (3) These five ESs will show concurrent degradation, and the degree 37 will be further intensified. (4) Future urban expansion will occupy more cropland and grassland 38 which will be the dominating reason for the intensified degradation of multiple ESs. 39 Conclusions We suggest that urban expansion through occupying a large amount of cropland and 40 grassland should be strictly controlled via urban land planning to alleviate the potential influences 41 of future urbanization on the concurrent loss of multiple ESs. 42


45
Ecosystem services (ESs) are the benefits that humans obtain from ecosystems (MEA, 2005). Urban 46 expansion refers to the process of transforming nonurban land into urban land under the expansion 47 of the urban scale in space (Bai et al., 2012). Urban expansion can promote regional socioeconomic 48 development and enhance the quality of human life, but it can also affect the structure and function 49 of ecosystems by removing vegetation and increasing the coverage of impervious surfaces, thus 50 leading to the concurrent degradation of ESs (Yang et al., 2020;Gong and Liu, 2021). Previous 51 studies have revealed that urban expansion has led to the concurrent loss of regulating services, 52 provisioning services, and supporting services ( world's land area, in which 38% of the world's population lives (MEA, 2005). With the rapid 57 expansion of urban land, the reduction in ESs is aggravated in drylands (Geist and Lambin, 2004). 58 In the future, drylands will still take place large-scale urbanization (Seto et al., 2012;Zhou et al., 59 2019). Therefore, the effective assessment of the potential influences of future urban expansion on 60 the concurrent loss of multiple ESs in drylands is of great importance for improving urban planning 61 and promoting regional sustainable development. in the northern part of Ningxia, China, and the results showed that urban expansion exerted negative 69 impacts on carbon storage, crop production, sand fixation, habitat quality, and nutrient retention. 70 Xie et al. (2018) assessed the impacts of urban expansion on water conservation, food production, 71 air quality regulation, habitat quality and carbon storage in Beijing, China, during 2013-2040 and 72 found that urban expansion will lead to the concurrent loss of multiple ESs in the future. However, 73 there are still some uncertainties in existing studies. This is mainly because the existing urban 74 expansion models have certain deficiencies in simulating future urban expansion. For example, the 75 urban expansion model used by Xie et al. (2018) is based on the traditional cellular automation (CA) 76 model. First, the traditional CA model assumes that the law of urban expansion remains unchanged 77 from the past to the future and uses static model parameters in the simulation. Second, errors from 78 data sources (such as field survey errors, digitization errors, and data conversion errors) will transmit 79 and accumulate constantly in the CA model, thereby increasing the uncertainty associated with 80 simulating future urban expansion (Zhang et al., 2015). 81 Improving the existing urban expansion model with the data assimilation method can lead to the 82 more accurate simulation of future urban expansion, thereby reducing the uncertainties of the 83 evaluation results of future urban expansion on the concurrent loss of multiple ESs. The data 84 assimilation method directly or indirectly integrates multiple sources and multiple resolution 85 observations by means of observation operators while weighing the uncertainties in the model and 86 the observations to minimize the error of the entire system (Sakov and Bertino, 2011;Li et al., 2020). 87 This method has been widely used in earth system sciences (Li et al., 2020 showed that urban land area increased from 151.29 km² in 1990 to 1,230.86 km² in 2017 in the 119 HBOY area, which caused the natural habitat quality to decline. With the development of the 120 regional economy and the relevant national policies (such as the plan for HBOY urban 121 agglomeration and the "One Belt and One Road" initiative), the HBOY will also undergo rapid 122 urban expansion in future, which will pose a huge threat to biodiversity and the surrounding 123 ecosystem (Song et al., 2020; National Development and Reform Commission, 2018). Therefore, 124 the HBOY urban agglomeration provides an appropriate study area for evaluating the potential 125 influences of future urban expansion on the concurrent loss of multiple ESs in drylands.

Data 135
In this study, the urban land data, land use/land cover (LULC) data, precipitation data, 136 socioeconomic data, and basic geographic information data were used.

Mapping ESs 152
According to the availability of data and the Millennium Ecosystem Assessment, we chose five key 153 ESs, including supporting services (natural habitat quality), provisioning services (food production), 154 regulating services (water retention and air quality regulation), and cultural services (landscape 155 aesthetic). is as follows: 161 where is the NHQ for grid j in natural habitat type i; is the suitability in natural habitat 163 type i; is the total threat for grid j in natural habitat type i; K is the half-saturation constant and 164 Z is the scaling parameter. The related parameters are obtained from Song et al. (2020). 165

Food production (FP) 166
Referring to Li

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

Landscape aesthetic (LA) 194
According to Cui et al. (2019), the visual quality index is used to quantify the LA, which represent 195 the aesthetic appeal of natural landscapes to tourists. The formula is as follows: 196 where is the total visual quality index score for grid j, with the range of 0-5; is the 198 terrain factor score for grid j; is the water bodies factor score for grid j; is the natural 199 vegetation factor score for grid j; ℎ is the urbanization factor score for grid j; and is the 200 natural landscape accessibility factor score for grid j. into urban grids is mainly affected by their suitability, the inheritance attributes of the LULC class, 209 neighborhood effects, and random perturbations in urban expansion. The probability for the 210 nonurban grid (x,y) with land-use type K transformed into the urban grid at time t can be expressed 211 as: 212

213
where ∑ −2 =1 × , , is the suitability of the nonurban grids (x,y) transformed into an urban grid 214 at time t; , , is the suitability factor i; is the weight of suitability factor i; , is the effects 215 of neighborhood grids; Wm-1 is the weights of neighborhood grids; , , is the inheritance attribute 216 of grid (x,y) at time t; is the ecological 217 constraints, which is a binary variable; ∏ , , =1 is the land planning policy constraints, which is also 218 a binary variable; and , is the random perturbation factor. 219 (2) EnKF model 220 The EnKF model can be divided into two stages: prediction and update. In the prediction stage, the 221 simulation of the process model is performed using the initial state set obtained by random sampling 222 to obtain the prediction set, and the Kalman gain is calculated by employing the prediction set. In 223 the update stage, the prediction set is updated by the observations and Kalman gain, and the updated 224 result is used as the analysis value set. The mean of the analysis value set is the posterior estimate 225 of the model and will continue to act as the new state set that participates in the next cycle (Sakov 226 and Bertino, 2011; Zhang et al., 2015). 227 The EnKF model achieves its goal of dynamic data assimilation by iterating the above two stages.

228
The main formula are as follows: 229 1) Initialization. The random variables Xi (i =1, 2..., N) were obtained by the Monte Carlo 230 method. In our study, the state variable is the urban land data in the initial year. 231 2) Prediction. The state variable at time k is used to estimate the state variable at time k+1. The 232 formula is as follows: 233 where +1 is the prediction set at time k+1; is the analysis value set at time k; and M() is the 235 state change relationship from time k to time k+1, which is generally a nonlinear operator. In our study, it is the LUSD-urban model.
represents the model error, which conforms to a Gaussian 237 distribution with an expected value of zero and a covariance of Q. 238

3) Update. The prediction set is corrected by integrating the observations and Kalman gain. 239
The formula is as follows: 240 where +1 is the analysis value set at time k+1; + is the Kalman gain at time k+1; +1 is 243 the observation set at time k+1; +1 is the observation operator at time k+1; and is the 244 observation error, which also conforms to a Gaussian distribution with an expected value of zero 245 and a covariance of R. ̅ +1 is the mean of the analysis value set at time k+1. 246 where +1 is the covariance of the prediction error at time k+1 and is the covariance of the 250 observation error at time k. 251 4) Judgment. If the end stage of assimilation has been reached, the model is terminated.

252
Otherwise, the model procedure returns to step 2). 253 (3) Combining the LUSD-urban model and the EnKF model 254 First, we used the factor data (i.e., distances to city centers, distances to highways, distances to rivers, 255 and distances to national roads), historical urban expansion data and the LUSD-urban model to build 256 the urban expansion model and obtain the initial parameters (Fig. 2). 257 states, namely, urban (cell state is 1) and nonurban (cell state is 0). According to Zhang et al. (2015), 262 we divided the study area into several grids, and each grid included a certain number of cells. We 263 calculated the urban expansion intensity of each grid as follows: 264 where , is the urban expansion intensity for grid (i,j) at time t; ∑ ( , = ) × is the 266 total number of urban cells in grid (i,j); and × is the total number of cells in grid (i,j).
where +1 is the analysis value set of weights at time k+1 and +1 is the prediction set of 277 weights set at time k+1. 278 where is the analysis value set of weights set at time k and is the error of weights, which 280 is represents Gaussian white noise. 281 Finally, we put the corrected results and parameters into the LUSD-urban model as the initial inputs 282 of the next cycle until all observations have been assimilated. After assimilation, the improved 283 LUSD-urban model can be used to simulate future urban expansion. 284 In our study, we put LULC in 2000, elevation, slope, distance to the city centers, distance to rivers, 285 distance to highways, distance to railways, and distance to national roads into the LUSD-urban which suggests that that the model is able to simulate future urban land expansion (Appendix C). 294 Appendix C Parameters used in the improved LUSD-urban model 296  t2. If ∆ , is more than zero, it indicates that ESs are enhanced. If ∆ , is negative, it indicates 321 that ESs are degraded and that urban expansion has a negative impact on ESs. 322 Then, Pearson's correlation coefficient was calculated to analyze whether there was concurrent loss 323 between the different ESs and quantify the size of the concurrent loss caused by urban expansion. 324 The calculation of the correlation coefficient is as follows: 325 where R is the correlation coefficient between ecosystem service X and Y; ( , ) is the 327 covariance between ecosystem service X and Y; ( ) is the variance of ecosystem service X; and 328 ( ) is the variance of ecosystem service Y. 329 When a pair of ESs revealed a trend of deterioration along with urban expansion (∆ , < 0), they 330 also had a significant and strong positive linear relationship (R >0.5, P < 0.01) both before and after 331 urban expansion, which indicates that urban expansion caused these two ESs to be concurrently lost. 332 In addition, the change in the correlation coefficient represents the change in the magnitude of the 114.38 km², and 114.07 km², respectively. These cities will account for 72.54% of the total area of 346 newly increased urban land throughout the whole region (Fig. 3B). 347

The influences of urban expansion on the concurrent loss of ESs from 1990 to 2017 353
Urban expansion from 1990 to 2017 led to a decline in the five ESs in the HBOY region. Among 354 them, the loss of FP was the largest, which was reduced from 132.87 t/km² to 131.77 t/km², for a 355 total decrease of 0.83%. The loss of LA was the smallest, which was reduced from 3.85 to 3.84, 356 decreasing by 0.26%. WR, AQR and NHQ decreased by 0.54%, 0.43% and 0.39%, respectively 357 (Fig. 4). Spatially, urban expansion in Baotou and Hohhot had the most serious impacts on the 358 evaluated ESs (Fig. 4). Urban expansion in Baotou from 1990 to 2017 led to reductions in NHQ, 359 FP, WR, AQR, and LA by 8.09%, 14.98%, 9.50%, 7.63%, and 5.17%, respectively (Table 1). Hohhot 360 decreased by 4.60%, 12.32%, 7.02%, 4.51% and 4.15%, respectively (Table 1). 361 2017 were greater than 0.5 and passed the significance level at 0.01 (Fig. 5a, b). Among them, the 365 correlation coefficients of NHQ and AQR were the largest in 1990 and 2017 and were both greater 366 than 0.8, which indicated that the concurrent loss of NHQ and AQR caused by urban expansion was 367 most severe in 1990-2017. In addition, the Pearson's correlation coefficients of the three pairs of 368 ESs showed an increasing trend, indicating that the intensity of the concurrent loss caused by urban 369 expansion was increasing. Among them, the intensity of the concurrent loss of WR and LA increased 370 the most, and the correlation coefficients increased from 0.57 to 0.62, representing an increase of 371 0.05. WR and FP, AQR, and NHQ increased by 0.03 and 0.01, respectively (Fig. 5a, b). 372    With the expansion of urban land in future, the five key ESs will continue to degrade. From 2017 to 393 2050, FP, WR, NHQ, AQR, and LA will be reduced by 1.34%, 0.81%, 0.65%, 0.60%, and 0.35%, 394 respectively, which will be greater that the reduction caused by urban expansion from 1990 to 2017. 395 Among them, the loss of FP will be the most serious and decreased from 131.77 t/km² to 130.01 396 t/km². The loss of LA was the smallest and decreased from 3.84 to 3.83 (Fig. 6). Spatially, the 397 influences of urban expansion on regional ESs in Baotou and Hohhot will still be the most serious 398 (Fig. 6). Among them, urban expansion in Hohhot from 2017 to 2050 will lead to reductions in FP, 399 WR, NHQ, AQR, and LA by 8.91%, 26.75%, 13.55%, 7.99%, and 6.56%, respectively (Table 2). 400 Those of Baotou will decrease by 11.52%, 20.65%, 12.04%, 8.94% and 5.19%, respectively (Table  401 2). 402 Urban expansion will lead to the concurrent loss of five pairs of ESs from 2017 to 2050, and the 403 intensity will become increasingly serious. Among them, NHQ and AQR, FP and WR, and WR and 404 LA will still show a concurrent loss. At the same time, FP and LA and WR and AQR will also reveal 405 concurrent losses in future. Specifically, the related ESs pairs will reveal a significant strong positive 406 correlation in 2017 and 2050 (Pearson's correlation coefficients will be greater than 0.5 and pass 407 the significance at 0.01) (Fig. 5b, c).The correlation between NHQ and AQR will still be the highest, 408 which indicates that the concurrent loss caused by urbanization will still be the most severe in the 409 future. In addition, with the advancement of urban land in the future, the Pearson's correlation 410 coefficients of the five pairs of ESs will also increase, indicating that urban expansion from 2017 to 411 2050 in the HBOY region will cause the concurrent loss of the five pairs of ESs, which will become 412 increasingly serious. Among them, the intensity of the concurrent loss between WR and AQR will 413 increase the most, and its correlation coefficient will increase from 0.54 to 0.59, representing an 414 increase of 0.05 (Fig. 5b, c). 415 Urban expansion will aggravate the loss of ESs and the degree of the concurrent loss of ESs in the 416 future. With the increase in urban areas, the average annual loss of the five evaluated ESs in 2017-417 2050 will be significantly higher than that in 1990-2017 (Figure 7). The most significant change in 418 the loss of ESs will be that of FP. Compared with 1990-2017, the average annual loss of FP in 2017-419 2050 will be 0.04%, representing an increase of 32.50%. In addition, compared with 1990-2017, the 420 pairs of ESs exhibiting concurrent loss will increase from three pairs to five pairs, and the correlation 421 coefficients will continue to show an increasing trend (Fig. 5).   was found that the Kappa coefficient of the improved LUSD-urban model was 0.64, which was also 455 higher than the 0.59 value output for the original model (Fig. 8). 456 Table 2    The assessment results using the improved LUSD-urban model and the ES models were closer to 465 the true value. Compared with the original LUSD-urban model, the results of combining the 466 improved LUSD-urban and the ES models showed that urban expansion reduced NHQ, FP, WR, 467 AQR, and LA by 11.99%, 50.89%, 26.77%, 17.94%, and 11.72%, respectively, which were closer 468 to the influences of urban expansion on ESs from 2010 to 2017 in the HBOY region (Table 3). In 469 addition, according to He et al. (2016), the relative error was also selected to further verify the 470 effectiveness of our method. The outcomes revealed that the relative errors of our method were -471 24.64%, 44.94%, 18.66%, 1.47%, and 8.02%, respectively, which were lower than those of the 472 original LUSD-urban model, and the simulation accuracy was improved by more than 20% on 473 average (Table 3). 474

Policy implication 481
Our study shows that regional urban expansion will ineluctably continue to cause the concurrent 482 loss of multiple ESs, especially for FP and WR. Urban expansion led to a 0.83% decrease in FP 483 from 1990 to 2017, with the highest impact reaching 14.98% in Baotou, followed by a 0.54% 484 decrease in WR, with the highest impact reaching 7.02% in Baotou. Urban expansion will lead to a 485 1.34% reduction in FP during 2017-2050, with the highest impact of 26.75% in Hohhot, followed 486 by a 0.81% decrease in WR, with the highest impact of 13.55% in Hohhot. In addition, both FP and 487 WR showed concurrent losses during these two periods, and the intensity will further expand in the 488 future. The large loss of FP will certainly bring great challenges to regional food security, and the 489 large decline in WR caused by urban expansion will certainly become the main bottleneck of urban 490 development in drylands. Therefore, the protection of FP and WR services should be emphasized in 491 future urban expansion in the HBOY region. 492 A large amount of cropland and grassland occupied by urban expansion is the dominating reason 493 for the concurrent loss of regional key ESs. We found that the loss of NHQ, FP, WR, AQR, and LA 494 caused by the occupation of cropland and grassland by urban expansion from 1990 to 2017 495 accounted for more than 70% of the total loss of each service (Table 4). This proportion will rise to 496 more than 75% by 2017-2050 (Table 4). In addition, urban expansion in the future will occupy more 497 cropland and grassland, which will be the dominating reason for the enhanced intensity of the 498 concurrent loss of the five ESs. We found that compared with the period from 1990 to 2017, the 499 amount of cropland and grassland occupied by urban expansion will increase by 47.38% from 2017 500 to 2050, which will lead to an average 1.44-fold increase in ES losses. 501 In the future, the local government should strictly control the new urban land occupying a large 502 amount of cropland and grassland through urban planning and management. For example, the 503 policies of "red line of arable farmland" and "ecological redline" should be strictly implemented in 504 the development of urbanization, or the protection of the high-quality cropland and grassland 505 ecosystem should be strengthened to alleviate the potential influences of future urbanization on the 506 concurrent loss of multiple ESs. In addition, future efforts should also aim to control the scale and 507 optimize the spatial pattern of cities to cut down the occupation of cropland and grassland ecosystem 508 by urban expansion. 509

Future perspectives 514
Our study proposed a new approach to effectively assess the potential impacts of future urban 515 expansion on the concurrent loss of multiple ESs in drylands. We combined the improved LUSD-516 urban model with ES models and used the HBOY urban agglomeration to verify the method 517 effectiveness. The results revealed that our method improved the accuracy of the assessment by 518 more than 20% on average, among which the accuracy of assessing the influences of urban 519 expansion on WR and AQR improved the most by 61.19% and 29.12%, respectively. 520 There are also some uncertainties in this study. First, the linear regression model was simply used 521 to determine the future urban population and urban land demand. However, the relationships 522 between economy, population, environment, policy and urban expansion are very complex. Second, 523 the quantification of ESs supply based on LULC data in this study was relatively simple, which may 524 lead to the underestimation of the impacts of urban expansion on the studied ESs. For example, the 525 influences of urban land on ESs due to the occupation of water bodies were ignored in this study. 526 However, this does not impact the dependability of the results of our study. This is because these 527 quantitative methods have been fully applied and verified in existing studies (Yang et al., 2015;528 Landuyt et al., 2016; Xie et al., 2018), and the occupation of water bodies by regional urban 529 expansion is considered to be relatively minimal, so the influences of urban expansion on the 530 concurrent loss of multiple ESs can be accurately depicted in our study. 531 In the future, we will use the mechanism model (such as system dynamics model) to predict the 532 future urban land demand based on the analysis of the driving mechanism of urban expansion. In 533 addition, we will further explore more effective and accurate methods of mapping ESs, thereby 534 reducing the uncertainties of the assessment results. 535

536
Our method can be used to effectively assess the potential influences of future urban expansion on 537 the concurrent loss of multiple ESs in drylands. First, the improved LUSD-urban model was 538 obviously superior to the original LUSD-urban model in the urban land simulation. The Kappa value 539 of the improved LUSD-urban model in simulating urban land in 2017 was 0.64, which was higher 540 than the 0.59 value of the original model. Second, compared with the original model, the influences 541 of regional urban expansion on multiple ESs evaluated by combining the improved LUSD-urban 542 model and ES models were closer to the actual value, with smaller relative errors. The accuracy of 543 the assessment was improved by more than 20% on average. Among them, the accuracy associated 544 with evaluating the influences of urban expansion on WR and AQR services improved the most at 545 61.19% and 29.12%, respectively. 546 Under the influence of urban expansion from 1990 to 2017, FP, WR, AQR, NHQ, and LA in HBOY 547 showed a downward trend, which decreased by 0.83%, 0.54%, 0.43%, 0.39%, and 0.26%, 548 respectively. During 2017-2050, urban expansion will accelerate the decline of these five key ESs, 549 and the average annual loss will reach 4.06‰, 2.47‰, 1.81‰, 1.97‰ and 1.05‰, which will be 550 higher than that of the 3.06‰, 2.01‰, 1.61‰, 1.45‰ and 0.99‰ losses that occurred during 1990-551 2017. At the same time, future urban expansion will lead to the concurrent loss of five pairs of ESs, 552 and the magnitude of the concurrent loss will be further intensified. The correlation coefficients 553 between them will increase from 0.52-0.86 to 0.57-0.87. Urban expansion in the future will occupy 554 more cropland and grassland, which will be the dominating reason for the concurrent loss of multiple 555 ESs. Compared with 1990-2017, the cropland and grassland occupied by urban expansion will 556 increase by 47.38% in 2017-2050, which will lead to an average increase of 1.44 times the loss of 557 ESs. 558 Therefore, we suggest that the future urban expansion encroaching on large amounts of cropland 559 and grassland should be strictly controlled in urban planning and management to reduce the risk of 560 the impact of urban expansion on multiple regional ESs, especially for cities undergoing rapid urban 561 expansion and serious impacts on ESs, such as Hohhot and Baotou. 562