Land use management based on multi-scenario allocation and trade-offs of ecosystem services in Wafangdian County, Liaoning Province, China

Developing effective methods to coordinate the trade-offs among ecosystem services (ES) is important for achieving inclusive growth and sustainable development, and has been the focus of scholars and ecosystem managers globally. Using remote sensing and geographic information system (GIS) data, our study examined Wafangdian County of Liaoning Province as a case study to reveal the spatiotemporal evolution of four ES (food supply [FS], net primary productivity [NPP], water yield [WY], and soil conservation [SC]) and changes among their interactions. Then, an ordered weighted averaging model was introduced to simulate the optimal scenario of ES allocation. Results showed that: (1) the spatial and temporal changes in ES were significant over 14 years. All ES presented an inverted U-shaped growth curve from 2000–2014. (2) Synergies were observed within provisioning services, and there were trade-offs between provisioning services and regulating services, as well as provisioning services and supporting services. (3) The optimal scenario for Wafangdian was scenario 5 (trade-off coefficient, 0.68). The allocation of FS, NPP, WY, and SC in scenario 5 were 0.187, 0.427, 0.131, and 0.063, respectively. Implementing each ES weight of optimal scenario in land use management contributed to achieving intercoordination of ES. We propose to coordinate land and sea management to restore natural habitats that were expanded into in the high ES area. It is our anticipation that this study could provide a scientific basis for optimizing the allocation of ES and improving land use structure of coastal zones in the future.

139 Where P i refers to the grain yield and animal husbandry at grid i. P sum represents the gross grain 140 yield and animal husbandry. NDVI i is the NDVI at grid i; NDVI sum is the sum of NDVI of 141 cropland or grassland in the study area. In addition, the yield of freshwater aquaculture and 142 marine aquaculture was assigned to inland waters and marine aquaculture areas, respectively.

143
Net primary productivity: NPP is an important index for evaluating vegetation 144 productivity and coverage and is an important component of the terrestrial carbon cycle (Field et 145 al., 1998). It was used to express the organic matter energy accumulated per unit area of 146 vegetation in a unit time; i.e., the residual amount of organic matter created by deducting oxygen 147 respiration from vegetation during photosynthesis (Ruimy et al., 1994). In our study, NPP was 148 modelled using the process-based Carnegie-Ames-Stanford approach (CASA) (Potter et al., 149 1993). Two important variables were used to account for NPP in the CASA model: the absorbed 150 photosynthetically active radiation (APAR) and the actual solar energy utilization efficiency of 151 plants (ε). 159 Where Y(x) is the annual WY of grid x (mm); AET(x) is the actual annual evaporation of grid x 160 (mm); P(x) is the annual precipitation of grid x (mm).

161
Soil conservation: Soil conservation was simulated using the revised universal soil loss 162 equation (RUSLE) model (Renard et al., 1997). Based on the replacement hypothesis of surface 163 cover, SC amount can be expressed as the difference between potential soil erosion and actual 164 soil erosion; i.e., the difference between the amount of soil erosion without any vegetation cover 165 management and soil and water conservation measures and the amount of soil erosion under 166 current vegetation cover management and soil and water conservation measures (Wu et al., 167 2017). The calculation formula used was as follows: 168 Where A c represents SC (t/hm 2 •a), A p is the amount of the potential soil erosion (t/hm 2 •a), and A r 169 is the amount of the actual soil erosion (t/hm 2 •a). The potential soil erosion and actual soil 170 erosion were quantified as follows: Where R is the rainfall erosion factor; K is the soil erosion index; LS represents the combined 172 slope length (L) and factor of slope length (S); C is the vegetation cover index; and P is the SC 173 measure index. Based on the correlation analysis, the trade-offs relationships among ES were clarified. 184 When a correlation coefficient among ES passed the significance test at 0.01 level and the 185 coefficient was negative, then there was deemed to be a significant trade-offs relationship. On 186 the contrary, if the correlation coefficient was positive, then there was a significant synergies 187 relationship . In our study, all the significant correlation coefficients were 188 classified into strong correlations (|r| ≥ 0.5), moderate correlations (0.3 ≤ |r| < 0.5), and weak 189 correlations (0.1 ≤ |r| < 0.3) according to the classification of the intensity of significant 190 correlation coefficients of Cohen (1992).

192
Ordered weighted averaging is a method for controlling factor weight combination. By 193 reordering the data of each index according to their attribute value and giving different order 194 weights and weighting aggregation values according to the order of each index, OWA reflects 195 the decision-making results when the decision maker ranks the importance of each index 196 differently (Yager, 1996;Zhang et al., 2017). Its model principle is as follows:

208
(2) The order weights w and trade-offs degree with different risk-degree coefficients were 209 calculated (Yager, 1993;Zhang et al., 2015) as follows (Table 1): Where w i is the order weight; Q RIM is a monotone increasing function; i is the ordinal number, 211 and n represents the number of ES; r is the independent variable; and α is the risk coefficient 212 (scenario), also known as the optimistic degree coefficient, which based on the risk perception of 213 decision-making caused by the difference in the numerical value of indicators and the difference 214 in subjective weights.

215
(3) Combining the OWA operator with the GIS platform, the OWA-based ES raster was 216 obtained.

217
(4) Using the Getis-Ord Gi* module in ArcGIS 10.2 to analysis the cold and hot spots of the 218 OWA-based ES raster under different scenarios. In our study, 99% of the hot spots were defined 219 as the highest ES area, 95% as a high ES area, and 0 as a medium ES area. On the contrary, 99% 220 of the cold spots were named as the lowest ES area, and 95% were low ES area. Finally, the 221 trade-off coefficient and the area of the highest and high ES areas of each scenario were 222 combined to determine the optimal trade-offs of ES in Wafangdian.

225
The maps in Figure 2 show the distribution of ES in 2000, 2007, and 2014. FS was highly 226 distributed in the central and western areas but had a smaller distribution in the east. The high 227 value areas were mainly in the vast central cropland and southwestern coastal waters, and the 228 overall trend decreased from west to east. The average FS yield in 2000 was 173.33 t/km 2 . In 229 comparison, the mean FS yield in 2007 was 250.57 t/km 2 , with a 44.56% increase compared with 230 2000. The mean FS production was 228.17 t/km 2 in 2014 and decreased by 8.94% compared 231 with 2007. In terms of temporal scale, FS yield showed a trend of rapid growth in the early stage 232 followed by a slow decline in the later stage. Among them, the FS capacity increased 233 significantly in the southwest coastal waters and the cropland along the lower reaches of the 234 Fuzhou River.

235
The total NPP of Wafangdian in 2000, 2007, and 2014 was 1.62 trillion g C, 1.85 trillion g 236 C, and 1.65 trillion g C, respectively (Table 2), showing a similar trend to that of FS in temporal 237 scale. However, the spatial distribution was different from that of FS. The high value areas were 238 mostly concentrated in the mountains with high vegetation coverage, such as Laomao Mountain 239 and Longtan Mountain in the east and Dabei Mountain in the middle and south. Low value areas 240 were mostly distributed in the southwest coastal areas, and in Fuzhou Town and the downtown 241 of Wafangdian.

242
The spatial pattern of WY was mainly influenced by precipitation and evapotranspiration. It

248
The average SC in 2000, 2007, and 2014 was 61.46 t/hm 2 , 66.28 t/km 2 , and 48.03 t/hm 2 , 249 respectively (Table 2), showing a moderate increase and then a sharp decrease through time. The 250 spatial distribution of SC was similar to that of NPP. Forest dominated the mountainous area in 251 the northeast, with high vegetation coverage; therefore, the SC had a high capacity. The terrain in 252 the southwest was flat, where towns and villages were concentrated. Owing to the high intensity 253 of land development and utilization, the SC capacity was low. On the whole, except for the 254 sustained growth of SC, the ES of the study area increased initially and then decreased from 255 2000 to 2014.
256 Spatiotemporal changes in ecosystem services with different land use 257 types 258 Our study focused on the changes in ES and their relationships among forest, grassland, and 259 cropland (three vegetation-covered land use types). The mean values of the four ES according to 260 each land use type were determined. To ensure a more intuitive comparison among the different 261 ES, min-max normalization was used to remove dimension and then Nightingale's rose diagram 262 was made. As shown in Figure 3, compared with cropland and grassland, forest had the largest 263 NPP and SC, while grassland had the largest FS. The WY of cropland was the highest, while SC 264 and NPP were the lowest, and FS was slightly lower than that of grassland. In terms of temporal . 271 Except for SC, the other three ES in cropland showed a significant inverted U-shaped trend 272 representing an initial increase that then decreased over the 14 years investigated (Fig. 3). 273 Overall, our results showed that NPP, FS, and WY decreased in forest, farmland, and grassland 274 from 2000 to 2014. SC in the three land use types increased slightly; however, this was not 275 significant (Fig.3).  Table 1. The OWA-based ES spatial distribution was different under each scenario. Compared 311 with the other six scenarios, the ES of scenario 1 presented low values. The spatial distribution of 312 the high and low values of scenarios 2-4 were more uniform than those of scenario 1. However, 313 the high values were more concentrated in the coastal mariculture areas, while the ES of 314 terrestrial areas were distributed evenly but were still dominated by low values. The trend in high 315 value clustering of ES under scenarios 5-7 became increasingly evident. In particular, there were 316 large sea-land differences under scenarios 6 and 7. The high values were mostly concentrated in 317 the eastern mountainous areas covered by forest, while the coastal waters were completely 318 covered by low values. Overall, from scenario 1-7, the spatial distribution of OWA-based ES 319 tended to be low-value clustering to uniform distribution to high-value clustering. 320 Scenarios 1 and 7 were two extreme scenarios (extremely optimistic and pessimistic, 321 respectively) in decision making, which were dominated by a single type of ES. These two 322 scenarios do not exist in reality and also do not meet the needs of diversification of ES. 323 Therefore, we excluded scenarios 1 and 7, and then the cold hot spot analysis module in ArcGIS 324 10.2 was used to determine the high (hot spot) and low (cold spot) value aggregation under the 325 other five scenarios (Fig. 7). For the spatial distribution of ES regionalization, scenario 2-6 326 changed from low value agglomeration dominant (scenario 2) to coexisting high-value and low 327 value agglomeration (scenario 6). Under scenario 5, the high ES areas radiated to the western 328 plains and to the coastal areas, including the highest ES observed. The smallest sea-land 329 differences in ES were observed under all scenarios.

330
The results of ES regionalization (Table 3) showed that the sum of the high ES areas under 331 scenario 5 covered an area of 2086.50 km 2 , accounting for 53.91% of the total area of 332 Wafangdian, far exceeding that of the other scenarios. Combining the trade-off coefficients of 333 each scenario in Table 1, scenario 5 had a trade-off coefficient of 0.68, which was the highest 334 under all feasible decision-making scenarios. Although the trade-off coefficient of scenario 4 335 was 1, the weights of services were equal under this scenario, which is an ideal decision-making 336 scenario for the complete equalization of ES. By comparing with the spatial distribution and area 337 of the high ES areas of each scenario and combining with the trade-off coefficients of each 338 scenario, it was found that scenario 5 had the largest high ES area (53.91% of the whole region) 339 and the highest ES trade-off coefficient (0.68). Consequently, the optimal trade-off scenario in 340 Wafangdian was the area defined under scenario 5. 341 Figure 8 shows that, in the optimal trade-off scenario, the high and highest ES areas of 342 Wafangdian covered all kinds of land use types and were distributed throughout the study area. 343 Cropland encompassed an area of 730.28 km 2 and constituted 35% of all of the high ES areas. 344 Forest and mariculture area were ranked second and third, respectively. Forest encompassed an 345 area of 542.49 km 2 and constitutes 26% of all of the high ES areas.
346 Discussion 347 Ecosystem services mapping and drivers of variation 348 In the present study, for the purpose of making better decisions on the optimization of 349 environmental and social services, methods based on remote sensing data were used for spatially 350 visualizing ES to enable stakeholders and policy-makers to compare with the actual land use 351 (Mina et al., 2017). Therefore, three ES use an evaluation model based on remote sensing or GIS 352 (SC based on RUSLE, WY based on InVEST, and NPP based on CASA). Maps were 353 comparable after standardization, which also facilitated the subsequent implementation of map 354 algebraic overlay and GIS spatial analysis. The maps produced by the models were the products 355 of a linear combination of multiple variables, such as evapotranspiration, soil composition, 356 vegetation cover index, and even statistical data of grain yield, which serve to supplement and 357 enrich land use data. The information provided by combinatory ability of the model was more 358 accurate and reliable than that obtained from the direct interpretation of land use (e.g., Eigenbrod 359 et al., 2010; Kandziora et al., 2013).

360
The map in Fig. 2 shows the spatial distribution of the four ES analyzed. The two 361 provisioning services (FS and WY) maintained the same spatial distribution pattern. This was 362 different from the results of Sun et al. (2016), whose results suggested that, where FS is high, the 363 water consumption of crops is also high, and thus, the higher the FS the lower the WY. These 364 results are reasonable for inland areas or research areas excluding offshore waters. However, the 365 large area used for aquaculture, and the highly important position of fisheries in the local 366 economy, introduce difficulties when attempting to classify coastal aquaculture areas into land 367 use types. The high production value and yield of mariculture areas result in a broader 368 distribution of FS in aquatic areas compared with that of terrestrial areas.  Table S1). These meteorological factors played a decisive role in the CASA, water yield module 377 of InVEST, and RUSLE models (Huston & Michael, 2012), which directly contributed to 378 increases in ES, such as NPP, WY, and SC, in 2007. In addition, our results clearly demonstrate 379 that the particular hydrothermal conditions in 2007 resulted in a higher crop yield in that year 380 compared with those of the same period over 14 years. Especially for Wafangdian, which is 381 located in the temperate monsoon climate zone, timely and sufficient precipitation and heat in 382 spring and summer (April-September) ensure the growth of crops and other plants. 383 Interaction among ecosystem services 384 The specific interactions between different ES were complex. The results of correlation 385 analysis among ES in different research areas are not consistent, especially regarding the 386 relationships among multiple ES in some areas where the interface between humans and land is 387 more prominent, which have not been clarified . In addition, the trade-offs 388 relationship has spatiotemporal heterogeneity, and the synergies relationship determined in one 389 region is often expressed as a trade-off in another region. For example, the study by Enfors et al. 390 (2008) at catchment scale showed that FS and SC are synergistic, while a study by Maes et al. 391 (2012) across the European continent showed a trade-off between these two ES. Therefore, the 392 trade-offs and synergies among ES derived from our study are not globally uniform, or even 393 nationally, and are also different from the results of previous studies.

394
Our study identified the synergies within provisioning services, the trade-offs between 395 provisioning services and regulating services, and the synergies between supporting services and 396 regulating services in Wafangdian. Previous studies have reported the synergies between WY 397 and regulating services, but trade-offs between WY and FS (Qiu & Turner, 2013;Sun et al., 398 2016). Conversely, our results revealed a more differentiated picture. One of the reasons for the 399 difference is as mentioned in the analysis of WY spatial distribution pattern in section 4.1; i.e., 400  Manuscript to be reviewed