The Weihe River, originating in Niaoshu Mountain in Weiyuan County, Gansu Province, has a basin area of 134,800 km² and an average annual runoff of 10.37 billion m³. The river enters Shaanxi Province from Baoji City, crosses the Guanzhong Plain and empties into the Yellow River in Tongguan County, Shaanxi Province, China. The average precipitation is approximately 550 mm, the altitude is approximately 240-3700 m, the regional topography is high in the west and south, and low in the east and north. Additionally, most of the region in the west, north and south is composed of mountainous areas, but the middle section is composed of the Guanzhong Plain, and there are a variety of landforms and rich vegetation types in the Weihe River Basin in the Guanzhong-Tianshui Economic Region (hereinafter referred to as “WRGT”) (Figure 1). In recent years, with the gradual promotion of western development and the implementation of the strategy of building a strong western province, the WRGT has become a very important region in the spatial economic layout of Shaanxi Province and in the country as a whole. However, the rapid economic development and population growth has resulted in an increase in the demand for construction land, and cropland continues to decline; as a result, the contradiction between people and land is prominent. The unreasonable land use layout and the agricultural production mode of extensive cultivation and monoculture have caused soil erosion to be more serious, and the ecosystems of the river basin and the environment changed to varying degrees (Shang, 2007; Zhan et al., 2011). Shaanxi has become one of the provinces with the most serious soil erosion in China, where 80% of the cropland and 70% of the population are distributed in the soil erosion areas, and the area of soil erosion accounts for 46.6% of the total land area of the province. Therefore, it is of great significance to the study of how water conservation changes with land use change and how to optimize the pattern of water conservation in this area for eco-environmental construction, resource development and land management.

The data used in this paper mainly include (1) the basic geographic information data of the WRGT, including the administrative county, river, road, elevation, and soil type data, which were mainly from the National Geomatics Center of China; (2) land use data for 2000, 2005, 2010, and 2015 and vegetation type data in 2010 and 2015 from the Data Center of Resources and Environmental Sciences of the Chinese Academy of Sciences (http://www.resdc.cn); (3) meteorological data, including precipitation, temperature, wind speed, relative humidity and evapotranspiration data, which were from eight meteorological stations at Changwu, Baoji, Xi’an, Yaoxian, Tongchuan, Wugong, Huashan and Tianshui and obtained from the National Meteorological Information Center; and (4) social and economic statistical data, mainly including population, GDP, etc., which were obtained from the Statistical Yearbook of Shaanxi Province and Statistical Yearbook of Gansu Province.

In this section, we introduce BBN theory and development, water conservation principles, and BBN construction methods for water conservation based on existing data, which lays a foundation for the further prediction of water conservation services.

3.1.1 BBNs

In 1988, Pearl (1988) proposed the concept of the BBN, which is a graphical network based on probabilistic reasoning. BBNs are a directed acyclic graph (DAG) that consists of nodes that represent variables and directed edges that connect them. Each node contains discrete states of variables, probability distributions corresponding to discrete states and conditional probability tables (Landuyt et al., 2016). The conditional probability table expresses the strength of the relationship between the parent node X and the child node Y. Each row in the table consists of the state combination of the parent node and the conditional probability P (y | x). In a network consisting of precipitation, vegetation type and evapotranspiration, the conditional probability table of evapotranspiration is shown in Table 1. Conditional probability tables can be obtained either from expert knowledge or from actual observation data. The probability distribution of the child node is determined by the probability distribution of the parent node and the conditional probability table. For a node without a parent node, the probability distribution is a prior probability distribution P(X). The joint probability of all variables in the BBNs can be obtained by multiplying the conditional probability distributions of all nodes (Formula 1). The joint probability characteristic enables the BBNs to effectively calculate the probability distribution of any node in the network. This result provides a theoretical basis for predicting the probability distribution of the water conservation node in this paper. This paper uses Netica software to construct a water conservation model based on BBNs; this software has powerful functions, an intuitive and smooth user interface, and allows users to easily draw and change the network.

$P({{X}_{1}},{{X}_{2}},\cdots ,{{X}_{n}})=\prod\limits_{i=1}^{n}{P({{X}_{i}}|parents({{X}_{i}}))}$ (1)

3.1.2 Quantification of water conservation

Currently, there are four main methods used to calculate water conservation: the soil water storage estimation method, the comprehensive water storage capacity method, the annual runoff method and the water balance method. The water balance method refers to the difference between the input and output of regional water. From the perspective of water balance, the difference between precipitation and forest evapotranspiration and other consumption is equal to the value of water conservation (Xiao et al., 2000). The method is the basis for studying the water conservation mechanism, as it fully reflects the distribution of regional precipitation and can more accurately quantify water conservation. Therefore, this study chose the water balance method to estimate water conservation in the WRGT; additionally, this method was used to choose the related variables as the nodes of the BBNs to construct the water conservation network model. Based on the data of vegetation types in the WRGT, the amount of water conservation in different land use types was calculated, and the amount of total water conservation was obtained after accumulation:

$Q=10\times \Sigma {{S}_{i}}({{P}_{i}}-{{E}_{i}}-{{R}_{i}})$ (2)

where Q is the total amount of water conservation, unit t/a; S_i is the area of land use type i, unit ha; P_i is the annual precipitation of land use type i, unit mm/a; E_i is the annual evapotranspiration of land use type i, unit mm/a; and R_i is the surface runoff of land use type i, unit mm/a.

3.1.3 The BBNs for water conservation

According to the water balance principle, the influencing factors related to water conservation were selected, and they included precipitation, vegetation type, land use in 2010, soil type, evapotranspiration and surface runoff. Water conservation is mainly influenced by

precipitation, evapotranspiration and surface runoff. Evapotranspiration is affected by precipitation and vegetation type (Cheng et al., 2014). Surface runoff is related to precipitation, land use and soil type (Cheng et al., 2014). Therefore, the Bayesian concept network for water conservation is constructed as shown in Figure 2. The raster layer corresponding to the node is discretized by ArcGIS software. According to the actual distribution of the variables in the study area, the raster layer of each variable is divided into four grades: the highest, high, medium and lowest. To facilitate writing, these four grades are called states 1, 2, 3, and 4, respectively. For example, “precipitation = 2” means that the precipitation is between 681 and 817 mm/a, and the discretization criteria are shown in Table 2. Among them, the evapotranspiration data came from 2010 MODIS data with a 250 m spatial resolution. The surface runoff data were calculated by the SWAT model, and the water conservation data for 2010 were calculated by the water balance method. By extracting the values of the raster layers into a single raster layer and inputting the attribute table into the BBN, the probability distribution of all nodes in the Bayesian network for water conservation can be obtained for 2010. On this basis, the conditional probability table of each node in the network is calculated by the conditional probability formula.

3.2.1 Land use scenario

Land use forecasting methods mainly include the Markov stochastic process model, the system dynamics forecasting model, the artificial neural network forecasting model, the CA model, the CA-Markov model, and the CLUE-S model (Chen and Wang, 2005). The CA-Markov model integrates the quantitative prediction ability of the Markov model and the cellular automata (CA) model to simulate the spatial change of complex systems. It can quantitatively simulate the spatio-temporal patterns of land use. Therefore, the CA-Markov model is used to simulate the land use pattern under three scenarios in 2050. (1) The conversion area between land use types in 2000 and 2005 was used as the element of the Markov state transition probability matrix. (2) According to the characteristics of the land use types of cropland, forest, grassland, water area, urban land and unused land, the slope, distance to road, distance to city center and distance to water area were selected as the driving factors. According to the characteristics of land conversion, certain driving factors were selected to calculate the potential distribution of land conversion types, such as the potential conversion of cultivated land into forest land. (3) Using 2005 as the baseline year, land use in 2010 was predicted, and a kappa precision test was conducted using the actual land use data for 2010; subsequently, the accuracy of the land use prediction was quantitatively evaluated. (4) Based on the land use maps of 2000 and 2005, the land use patterns under different scenarios were simulated and predicted for 2050.

Previous studies have shown that the main driving forces of short-term land use change are social factors such as national policy, economic development and population growth (Gan et al., 2004). Combined with the data of the 2015 Statistical Yearbook for the study area, the scenario was constructed using the key factors of economic development and population growth. The specific indicators are shown in Table 3. Based on the key factor index values of each district and county under three scenarios, the spatial distribution of the index was realized by the interpolation method, and the potential distribution of land use type transformation under specific scenarios was generated; thus, the land use in 2050 was predicted. The planning scenario reflects the evolution of land use change according to its own development process, i.e., the future land use change naturally evolved according to the current economic development trend. The protection scenario reflects the overall land use planning of Shaanxi Province (2006-2020) to coordinate land use and ecological construction and to adjust land to local conditions. The task of improving the terrestrial eco-environment reflects the situation of restricting the speed of urban expansion, slowing economic development and protecting the eco-environment. The policy of urbanization reflects changes in land use types when the economy and the urbanization process are accelerating.

3.2.2 Water conservation prediction

The Bayesian network has a strong ability to address uncertainty. It can express the correlation between various informational elements based on conditional probability and the model can learn under conditions of limited, incomplete and uncertain information. Therefore, this paper combines the CA-Markov model with the Bayesian network model to construct a water conservation service forecasting model for 2050 by using the water conservation service network in 2010 and the land use types under different scenarios in 2050 (Figure 3), and the water conservation service network in 2010 is used as the first model. With the help of Netica and ArcGIS, the conditional probability table and the probability distribution of each node in the water conservation network in 2010 were obtained. The water conservation network under the scenario analysis in 2005 was the second model where the uncertainty of land use was introduced. The uncertainty associated with the distribution of each node in 2010 was immediately integrated into the conditional probability. The probability distribution of water conservation under the different scenarios in 2050 was obtained by Formula (1).

After the construction of the Bayesian network model of water conservation in 2050, the key state subset of key variables was determined by a conditional probability table and sensitivity analysis and was visualized by ArcGIS software. Finally, the key state subset of key variables corresponding to different states of the water conservation node was plotted. The determination of the critical state subset of key variables mainly included two processes:

(1) Determine the key state of variables.

Based on the BBNs, the probability of each node state and the joint probability of two node state combinations were derived, and the conditional probability was calculated between each influencing factor and water conservation. The state with the maximum conditional probability of the influencing factor corresponding to different states of the water conservation node was selected as the key state of the variable. When the state of the water conservation node was in the “highest”, “high”, “medium” and “low” states, the corresponding variable states were called “I”, “II”, “III”, and “IV”, respectively.

(2) Identify key variables.

To evaluate the relative importance of the Bayesian network nodes, Netica provides a sensitivity analysis (Landuyt et al., 2016), which assesses whether the water conservation nodes are sensitive to changes in the other impact factor nodes. Sensitivity analysis can be calculated by variance reduction (VR). The VR represents the difference between the variance of the variable ES and the variance of the variable ES under the condition of known variable I. The formula is as follows:

$VR=V(ES)-V(ES|I)$ $=\sum\limits_{s}{p(s)\times {{(s-E[ES])}^{\text{2}}}-\sum\limits_{s}{p(s|I)\times {{(s-E[ES|I])}^{2}}}}$ (3)

where s is the state of the output variables; ES is the output variables; I is the other variables in the network; P(s) is the probability value when the state of ES is s; P(s|I) is the conditional probability between ES and I; E[ES] is the mathematical expectation of ES, and E[ES|I] represents the mathematical expectation of ES when I happens. As the value of VR increases, the sensitivity of the water conservation nodes to the input nodes becomes stronger.

After determining the key state subset of the key variables of water conservation under the three scenarios, the scenario with the largest ratio of subset I area to the total area of the study area was selected as the optimal scenario for the water conservation ecosystem service. The key state subset I of the key variable under the scenario was identified as the optimal state subset. Based on the combination of the optimal state subset and the latest land use planning policy, the optimization zone was divided.

The BBN for water conservation in 2010 is shown in Figure 4. The probabilities of the highest, high, medium and low water conservation values are 0.151, 0.607, 0.206 and 0.037, respectively. The water balance formula was used to calculate the water conservation in2015. Two hundred data points were collected randomly in the study area. The data of precipitation, vegetation type, land use and soil type in 2015 were input into the network where all the data were located. The probability distribution of water conservation prediction was then calculated. The accuracy evaluation of the BBN for water conservation was conducted using an error matrix (Congalton, 1991). The error matrix between the predicted water conservation and the actual calculated water conservation was 75.5%, as shown in Table 4. The results show that the model has good recognition and strong reliability for the prediction of future water conservation.

When comparing the land use forecast map and the real land use distribution map for 2010, the kappa coefficient was 0.884, which showed that the land use simulation had high accuracy and credibility and could be used for future land use change predictions. IDIRSI software was used to simulate the three land use scenarios. The predicted land use change and distribution types are shown in Table 5 and Figure 5, respectively. Under the protection scenario, the reduced ratio of cropland area was 13.22%, the increased ratio of forest was 18.12%, the reduced ratio of grassland was 0.73%, and the reduced ratio of water area was 5.08%. The area of urban land increased by 0.38%, while that of unused land decreased by 0.92%. In the planning scenario, the area of cropland and water area decreased, and the area of forest increased. In the development scenario, the cropland was reduced by 11.52%, and the area of urban land increased by 5.04%. Compared with the three scenarios, the overall trend was the reduction of cropland and the increase of urban land and forest, which was in line with the background conditions of economic development in the study area. Under the protection scenario, the area of cropland decreased most, while the area of urban land increased most under the development scenario. Under the planning scenario, the areas of cropland, water area, grassland decreased, and the area of urban land increased. Although the forest area increased, the increase range was between that of the protection scenario and that of the development scenario, and the degree of urban expansion was higher than that of the protection scenario. Therefore, the overall eco-environmental quality of the planning scenario was in the middle of the three scenarios, and the eco-environmental quality under the development scenario was the worst, while it was the best under the protection scenario.

The land use under the three scenarios was used as the input data for the BBNs for water conservation, and the water conservation under each scenario in 2050 was obtained. The probabilities of highest, high, medium and lowest water conservation values under the protection scenario were 0.184, 0.574, 0.206 and 0.036, respectively; under the planning scenario, the values were 0.171, 0.580, 0.205 and 0.045, respectively; and under the development scenario, the values were 0.171, 0.565, 0.204 and 0.060, respectively. The probability of the water conservation state distribution under the three scenarios was similar, i.e., the probability of the water conservation value distribution was the highest in 917-1103 t/ha·a, and the probability of being lower than 163 t/ha·a was the lowest. At the same time, there were also differences, e.g., the probability of high water conservation value under the protection scenario was the largest, followed by the planning scenario and the development scenario, which indicated that land use under the protection scenario was more conducive to enhancing water conservation.

Sensitivity analysis of the BBNs for water conservation in 2010 was conducted, as shown in Table 6. The VR of precipitation was the highest, indicating that precipitation had the greatest impact on the final water conservation, followed by evapotranspiration and land use. The VR values of the other variables were less than 0.1, indicating that water conservation had a weak influence on the changes of these variables. Therefore, precipitation, evapotranspiration and land use were selected as the key variables used to determine the state of water conservation under different scenarios, i.e., the variables that have the greatest impact on water conservation. The conditional probability table between each node of the BBN and the water conservation node in different scenarios in 2050 was derived from Netica software. The results show the combination of various variable states and corresponding conditional probabilities in the corresponding position when the water conservation ecosystem service in a certain geographical location is in a different state. The key state subset of the three scenarios was the same: the key state subset of the key variables I = {precipitation = 1, evapotranspiration = 2, land use = 2}, the key state subset of the key variables II = {precipitation = 3, evapotranspiration = 4, land use = 1}, the key state subset of the key variables III = {precipitation = 4, evapotranspiration = 3, land use = 1}, and the key state subset of the key variables IV = {precipitation = 3, evapotranspiration = 1, land use = 5}.

The key state subset of the variables mentioned above is displayed spatially in Figure 6, and subset I is the set of variables with the highest probability of the highest water conservation capacity. The corresponding position of the subset has the characteristics of large annual average rainfall, large annual evapotranspiration and high vegetation coverage. It is mainly distributed in the south of Maiji District of Tianshui City, the south of Weibin District of Baoji City and the southwest of Longxian County, Linyou county, Fengxiang County, Qianyang County, the northeast of Xunyi County of Xianyang City, and the northwest of Yongshou County. The junction of Chunhua, Liquan and Jingyang counties and the western part of the Tongchuan Yaozhou District comprised subset II, which included a set of variables with a high probability of high water conservation capacity, corresponding to the areas with low annual average rainfall and low annual average evapotranspiration. The main land use type was arable land, which was mainly distributed in the northwest of Guanzhong Plain and Beishan Mountain, including the central part of Fengxiang County in Baoji City, Chunhua County in Xianyang City, Binxian County, Qianxian County, northern Sanyuan County, central Jingyang County and Liquan County, and the southeast of Yaozhou District in Tongchuan City. Subset III is the highest probabilistic variable set with medium water conservation capacity, which corresponds to areas with low annual average rainfall, low annual average evapotranspiration and high vegetation coverage. It is mainly distributed in northeastern Longxian, southern Fengxiang, Qishan, central Fufeng, Qianxian and southern Liquan counties. Subset IV is the most likely variable set with the lowest water conservation capacity. The corresponding region has lower annual average rainfall and large annual average evapotranspiration. The land use types are mainly found in Jintai District of Baoji City, as well as Fengxiang, Qianxian, Liquan and Sanyuan counties and Yaozhou and Wangyi districts of Tongchuan City. Subset IV is the highest variable set with the lowest probability of water conservation capacity, which corresponds to the lower annual average rainfall and the large annual average evapotranspiration. The land use types are mainly cities, which are primarily distributed in the districts and counties of the main urban areas as Jintai District of Baoji City, and Fengxiang, Qianxian, Liquan and Sanyuan counties, Yaozhou District of Tongchuan City as well as Wangyi District.

In the protection scenario, the area of subset I accounted for the highest proportion of the total area, i.e., 5.55%; in the planning scenario, the value was 5.27%; and in the development scenario, the proportion of the subset was 5.34%. Therefore, to construct the eco-environment of water conservation in the Guantian section of the Weihe River Basin, we should refer to the protection scenario and combine the land use orientation formulated by the General Land Use Planning of Shaanxi Province (2006-2020), that is, we should comprehensively renovate the Qianshan and Longshan hilly areas in northern Baoji City, effectively control soil erosion and nurture and protect forest land. Forest resources should be combined with agriculture, forestry and animal husbandry; eco-environmental construction should be strengthened in the hilly and mountainous areas in the north of Xianyang, the gullied area in the middle of tableland and the mainstreams of the Jinghe and Weihe rivers. Additionally, natural forest resources should be protected, protective forests should be built, silting-up dam harnessing projects should be popularized, ecology should be restored, and harmonious development between man and nature should be promoted. Finally, Yaozhou in Tongchuan should continue to develop. We will promote the comprehensive development and demonstration of the mountainous areas in Yaozhou and the demonstration counties (districts) of “changing slopes into terraces” to speed up the pace of land ecological construction. Therefore, to optimize the water conservation pattern, the south of Maiji District in Tianshui City, the southwest of Longxian County in Baoji City and the south of Weibin District, the northeast of Xunyi County and the northwest of Yongshou County in Xianyang City, and the west of Yaozhou District in Tongchuan City should be designated as the key optimization areas, and the construction of greening projects should be strengthened, as shown in Figure 7.

Land use change has an important impact on water conservation (Fan et al., 2017). This paper uses the CA-Markov model to simulate land use change. Using 2005 as the baseline year, the land use in 2010 was forecasted and compared with the actual land use in 2010. The kappa coefficient was 0.884, which indicated the high accuracy and feasibility of the land use simulation. Scenario analysis is an effective tool for exploring the risks of existing land use patterns and policy choices (Nelson et al., 2009; Guo, 2011). This paper considers the effects of natural and human factors on land use change. The slope, distance to the road, distance to the city center, and distance to the water area were selected as the driving forces of land use type transformation, and economic development and population growth were the limiting factors used to predict the future three land use scenarios. Using the CA-Markov model and social economic indicators to set up the land use scenarios, rather than simply setting node states as scenarios in the network, as was done previous Bayesian network research (Xue et al., 2017), the reliability of the scenario design was improved. However, the uncertainties of some driving factors affecting land use change and the new policies made in the future will lead to the repeated changes of land use types. How to overcome these uncertainties is a difficult problem to solve in land use change prediction in the future. In this paper, the land use forecasting model is coupled with the BBN, and the uncertainty of land use forecasting is transmitted to the water conservation node through a conditional probability table. Compared with other non-statistical models, the uncertainties of land use distribution under each scenario are considered to improve the accuracy of predicting the future distribution of water conservation.

Generally, the uncertainty of a BBN includes structural uncertainty, input data uncertainty and parameter uncertainty (Chen and Pollino, 2012). According to the water balance principle and the literature, this paper combined water conservation with BBNs to clearly and transparently show the process relationship of water conservation. Although the impact of human activities, such as economic development and population growth, were considered in combination with land use forecasting models, other factors, such as temperature, slope and reserve area policy, also have impacts on water conservation. Further research should consider adding more impact factors into the network to reduce the structural uncertainty caused by an incomplete understanding of water conservation processes, as this would make the results more practical. The input data in this paper are from the measured or calculated raster layer. There are no missing data, and the uncertainty is very low. The uncertainty of model parameters is caused by the state classification of the Bayesian network nodes. In this paper, the natural classification method of ArcGIS was used to classify the continuous variables in the network, and the literature and expert knowledge were used to classify the discrete variables (Li et al., 2019). In the future, more accurate methods could be used to classify node variables, such as dynamic discrete classification (Xue et al., 2017). In this paper, the overall accuracy of the BBN of water conservation was 0.755, which showed that the model was reliable and could be used to predict the distribution of water conservation in the future. However, although this study validated the model, it did not further calibrate the model. In the next study, the model should be repeatedly validated and calibrated with the above uncertainty analysis to continuously improve the model.

This study provided an open Bayesian network framework for water conservation, which was flexible in structure. In the future, it can also construct the relationship network between water conservation and other ecosystem services, including ecosystem functions, human activities, natural factors and policies, and it can be used to explore the trade-offs and synergies among ecosystem services under different scenarios and analyze the effects of driving factors on trade-offs and synergies among ecosystem services.

Currently, the research on the spatial pattern optimization of water conservation mainly focuses on establishing scenarios around land use pattern changes, forecasting the changes of water conservation under various scenarios, and providing optimization suggestions. However, land use is not the only factor affecting ecosystem water conservation, and the degree of impact of each factor on water conservation is also different. Based on BBNs, with the aid of GIS, this paper proposed a key variable key state subset method that combined scenario analysis to optimize water conservation. The results provide decision support for the delineation of water conservation eco-environmental protection zones in the study area and have certain reference value for the spatial pattern optimization of other ecosystem services.