Spatial-temporal pattern study on water conservation function using the SWAT model Open Access

Water resources are a fundamental resource for the sustainable development of agriculture, society, and economy in any region (Cao et al. 2020; Immerzeel et al. 2020). The water conservation function (WCF) is an essential function of an ecosystem to maintain water resources, including natural processes such as atmosphere, water, vegetation, and soil (Kurihara et al. 2018). Manifestations of the WCF are mainly the flood retention and runoff regulation of an ecosystem. Flood retention means that the ecosystem intercepts part of the rainwater through the forest canopy, litter layer, and soil layer of the vegetation, reducing the direct fall of precipitation on the surface, thereby significantly reducing the probability of flooding (Zhou et al. 2010; Návar 2017). Runoff regulation means that the ecosystem stores rainfall in the rainy season and replenishes rivers in the dry season, thereby stabilizing river flow in the dry season. The role of WCF in flood interception and runoff regulation are different manifestations of WCF at different time scales. The evaluation of WCF is usually represented in terms of the water conservation amount (WCA). The calculation methods of WCA are mainly the canopy interception method (Davies-Barnard et al. 2014), comprehensive water storage method (Liu et al. 2016), and water balance method (Scordo et al. 2018).The water balance method is the basis for studying the mechanism of water conservation, and the calculation results are more credible than other methods. It is currently the most commonly used and most extensive method in the quantitative study of WCFs. The water balance method regards the ecosystem as complex, and from the perspective of water balance, the difference between precipitation and evapotranspiration and other consumption is considered the WCA (Ding et al. 2017). In the current research, the water balance method and the ecohydrological model are combined to carry out quantitative analyses on WCF at the basin scale.

Affected by topography, climate, soil, vegetation, and other factors at different time and space scales, the WCF has a significant time scale effect and regional differentiation effect (Wang et al. 2013). Different regions and different time and space scales may contain completely different manifestations of WCFs. Therefore, a multi-time scale study of the spatiotemporal changes of the WCF can comprehensively evaluate an ecosystem's WCF. Most current studies use the average runoff over many years to calibrate model parameters (Cong et al. 2020; Xu et al. 2020). However, these studies can not reflect the inter-annual changes in the WCF, let alone the intra-year changes. Spatial stratified heterogeneity is one of the fundamental characteristics of geographic phenomena. At present, most studies on the factors affecting the spatial differentiation of WCFs consider the correlation or principal component analysis between WCA and climate, human activities, and other factors, and there is a lack of quantitative analysis of the influencing factors (Luo et al. 2019).

As an indispensable tool for understanding the interaction between hydrological processes and the ecological environment, the SWAT eco-hydrological model has been widely used in ecohydrological research at the basin scale (Yang et al. 2018). The SWAT model can accurately simulate the hydrological cycle's physical processes, such as interception, surface runoff, evapotranspiration, infiltration, and underground runoff from yearly, monthly, and daily time scales through the division of hydrological response units and sub-basins. Simultaneously, it can reflect the unevenness of precipitation, evapotranspiration, and spatial distribution of the underlying surface. Combining the SWAT model with the water balance method can provide a new way for the multi-time scale study of spatiotemporal WCF changes. GeogDetector is a new spatial analysis model that detects the spatial differentiation of geographic elements and reveals the influencing factors (Wang & Hu 2012). Geographic detectors can not only quantitatively analyze the weight of each factor but also analyze the interaction between multiple factors.

On April 1, 2017, China's Central Committee and State Council decided to establish the Xiong'an New Area to build a green, ecological, and livable city. The upstream of Xiong'an New Area is the most crucial WCF area and ecological environment support area of the city (Yang et al. 2018). The water conservation capacity of the ecosystem in the upstream basin dramatically affects the water supply and flood control safety of Xiong'an New Area. In this study, we take the upstream of Xiong'an New Area as the research area. Based on the basin SWAT model's construction, we combined the SWAT model with the water balance method to calculate the WCA at annual and monthly scales and analyze its temporal and spatial changes. Simultaneously, multiple linear regression and geographic detector models were used to study the influencing factors of WCA temporal changes and the influencing factors of the WCA spatial pattern. This study provides a theoretical basis for decision-makers to formulate and adjust the upstream water resources management policy and ecological environment protection policy in the upstream of Xiong'an New Area.

The study area is located in northern China's mountainous area (38° 46′ 31″ ∼ 39° 40′ 8″ N, 113° 40′ 4″ ∼ 115° 10′1″ E), and includes three basins in the upstream of Xiong'an New Area – Zijingguan (ZJG), Zhongtangmei (ZJG) and Fuping (FP) basins – with a total area of 7,310.27 km² (Figure 1). The area's altitude ranges from 149 to 2,611 m, gradually decreasing from northwest to southeast, and the main landform types are mountains and hills. The study area is located in a temperate monsoon climate zone, with hot and rainy summers, cold and dry winters, with an average annual temperature ranging from 7.4 to 12.8 °C and precipitation between 550 and 790 mm (Cui et al. 2019). Affected by the monsoon climate, the annual precipitation in the study area is uneven. Nearly 80% of the precipitation is concentrated in the rainy season from June to September (Zhao et al. 2012).

The data used in this study include meteorological data, such as precipitation and temperature, surface runoff data, a digital elevation model (DEM), land use, soil type, and other spatial data. The resolution, time, and sources of socio-economic data such as gross domestic product (GDP) (10⁴ yuan/km²) and population density (persons/km²) are shown in Table 1. Soil hydrological properties are calculated based on the soil water characteristic software SPAW (Soil Plant Atmosphere Water) developed by the United States Department of Agriculture. The locations of four meteorological stations and the 42 evenly distributed precipitation stations are shown in Figure 1.

The SWAT model divides the basin into multiple sub-basins according to the river network and water system, completes topological connections based on paths such as rivers, and then divides the sub-basins into smaller hydrological response units based on features such as surface cover, soil, and slope (Moriasi et al. 2007). The SWAT model used the Hydrological Response Unit (HRU) as the minimum simulation unit to simulate hydrological cycle processes. Previous research has shown that setting the area threshold of land use, soil, and slope to 0 can make the generated HRU reflect the actual amount and spatial distribution of land use (Kan et al. 2017).

In the SWAT model construction of this study, the surface runoff was calculated using the soil conservation service (SCS) runoff curve model, the soil water content simulation was made using the dynamic water storage reservoir model, the reference evapotranspiration calculation used the Penman-Monteith formula, and the groundwater runoff was calculated using the basal flow regression coefficient method (Han et al. 2012).

SWAT-CUP is a calibration software specially developed by the Swiss Federal Institute of Water Quality Science and Technology for the SWAT model (Khatun et al. 2018). The software realizes automatic parameter calibration and can give the optimal parameter value and the optimal parameter value range. The SUFI-2 algorithm in SWAT-CUP software is a common method for uncertainty analysis of hydrological models (Cao et al. 2018). The SUFI-2 algorithm takes into account uncertain sources such as model structure, monitoring data, and model parameters. Therefore, this paper intends to select the SUFI-2 method for watershed parameter calibration. The simulation results were evaluated using Nash-Sutcliffe efficiency (NSE) (Nash & Sutcliffe 1970), correlation coefficient (Zuo et al. 2016), and relative bias (PBIAS) (van Griensven et al. 2006). The model parameters and simulation results can be found in Wang et al. (2021). The results of model calibration and verification are shown in Table 2.

In this study, we first used ArcGIS 10.1 to perform spatial overlay analysis on the WCA layer and impact factor layer in 2010, and reclassified the factors to obtain type variables. The value of the dependent variable Y is a numerical quantity, representing the area attribute of each WCA patch. The independent variables X1, X2, …, X8 are type quantity, which are respectively the categorical attribute values of precipitation, potential evapotranspiration, altitude, slope, soil type, land use, GDP, and population density. We bring the X and Y data into the geographic detector model to get the q value result, and finally get the maximum q value by adjusting the classification criteria of the independent variable X. The value range of q is [0–1]. The larger the value, the stronger the explanatory power of this factor to the Y value (Todorova et al. 2016).

With the aid of the interactive detector module of the geographic detector, we evaluate whether the two impact factors will increase or decrease the explanatory power of the independent variable X to the dependent variable Y when they work together. The larger the value, the more it means that the two impact factors have a more significant effect on Y than a single impact factor.

It can be seen from Figure 2 that the inter-annual variations of WCA in the ZJG, ZTM, and FP basins were considerable. The multi-year average value of WCA was largest in the ZJG basin, which was 47.59 mm, and the smallest was in the FP basin, which was 4.55 mm. The maximum WCA in ZJG, ZTM, and ZTM basins all appeared in 2011, and the values were 47.59 mm, 103.88 mm, and 49.26 mm, respectively. At this time, most regions of the study area were positive regions of WCA. The high positive value regions of WCA mainly appeared in the ZJG basin and the middle and lower reaches of the ZTM basin; the moderate positive value regions of WCA mainly appeared in the southwest of the FP basin; the low positive value regions of WCA mainly appeared in the upstream area of the ZTM basin (Figure 3(e)). The minimum values of WCA in the ZJG and FP basins were both negative, and the values were −8.4 mm in 2014 and −55.65 mm in 2017, respectively. In 2014, most of the ZJG basin was a low positive value region of WCA. In the upstream of the ZJG basin, sub-basins 5, 8, and 10 showed a moderate negative value region of WCA, while sub-basin 7 showed a high negative value region of WCA (Figure 3(h)). In 2017, the FP basins were all areas with negative WCA, and the middle and lower reaches of the FP basin were dominated by a moderate and high negative value region of WCA (Figure 3(k)). The minimum value of WCA in the ZTM basin was positive, which was 1.86 mm in 2009. At this time, the upstream of the ZTM basin was mainly a moderate positive value region of WCA, while the upstream of the ZTM basin was mainly a low and moderate negative value region of WCA. Furthermore, a low positive value region of WCA showed in sub-basins 2, 3, 4, 7, 12, and 14 (Figure 3(c)).

It can be seen from Figure 4 that there were apparent differences in WCA in different months in each basin, which was not only reflected in the magnitude but also positive and negative changes. A negative WCA indicated that the WCA in the month was at a loss. The variability of WCA in the ZJG and ZTM basins was largest in July, and the variability of WCA in the FP basin was largest in August. The WCA in the rainy season (June-September) of each basin was positive, while the WCA in the non-rainy season (October-May) was negative.

The WCA of the ZJG, ZTM, and FP basins all reached their maximum in July, with values of 32.12 mm, 28.87 mm, and 37.18 mm, respectively (Figure 4). At this time, the study areas were all positive value regions of WCA, and the high positive value regions of WCA mainly appeared in the lower reaches of each basin. The moderate positive value regions of WCA mainly appeared in the upper reaches of the ZJG basin, the middle reaches of the ZTM basin, and the middle and upper reaches of the FP basin. The low positive value regions of WCA mainly appeared in the upper reaches of the ZTM basin (Figure 5(g)). The lowest values of WCA in ZJG, ZTM, and FP basins were all negative values. The lowest WCA values in the ZJG basin and FP basin occurred in November and October, respectively, with values of −8.61 mm and −17.89 mm, respectively. In November, the ZJG basin was dominated by moderate and low negative value regions of WCA, and the high negative value regions of WCA appeared in sub-basins 6 and 7 (Figure 5(k)). In October, the FP basin was mainly a high negative value regions of WCA, and only the upstream sub-basins 1 and 2 were moderate negative value regions of WCA (Figure 5(j)). The WCA in the ZTM basin was the lowest in May, with a value of −8.56 mm. At this time, the ZTM basin was mainly a moderate and low negative value region of WCA (Figure 5(e)).

We used the statistical package SPSS 22.0 to perform multiple linear regression to obtain the relationship and significance between the inter-annual variation of WCA and precipitation, potential evapotranspiration, and population density.

It can be seen from the R² results in Table 3 that the fitting equation can reflect the relevant factor information very well, indicating that there was an excellent linear relationship between WCA and the influencing factors. The WCA in the ZJG and ZTM basins showed a significant positive correlation with precipitation (p < 0.05), while the WCA in the FP basin showed a significant negative correlation with potential evapotranspiration (p < 0.05).

It can be seen from R² in Table 4 that there was an excellent linear relationship between WCA and influencing factors, and the R² of the monthly scale was greater than the R² of the annual scale in each basin. The WCA in ZJG and FP basins was significantly positively correlated with precipitation (P < 0.05). There was a significant positive correlation between WCA and precipitation in the ZTM basin (P < 0.05) and a significant negative correlation between WCA and potential evapotranspiration (P < 0.05).

The factor detector analysis results showed (Table 5) that the three most prominent factors affecting the spatial distribution of WCA in the ZJG basin were: potential evapotranspiration (q = 0.813) > precipitation (q = 0.783) > soil type (q = 0.306). The three most prominent factors affecting the spatial distribution of WCA in the ZTM basin were: precipitation (q = 0.891)> potential evapotranspiration (q = 0.817) > altitude (q = 0.506). Additionally, soil type also has a more significant impact on the spatial distribution of the ZTM basin, with a q value of 0.461. The three most prominent factors affecting the spatial distribution of WCA in the FP basin were: precipitation (q = 0.730) > potential evapotranspiration (q = 0.542) > altitude (q = 0.191).

The analysis results of the interaction detector show that the top three interaction factors of explanatory power in the ZJG basin were PRE ∩ PET (q = 0.953), PET ∩ Soil type (p = 0.932), PRE ∩ Soil type (q = 0.853) (Table 6); the top three interaction factors of explanatory power in the ZTM basin were PRE ∩ PET (q = 0.987)>PET ∩ Soil type (p = 0.944)>PRE ∩ Soil type (q = 0.935) (Table 7); the top three interaction factors of explanatory power in the FP basin were PRE ∩ PET (q = 0.993)>PRE ∩ Soil type (p = 0.860)>PRE ∩ Altitude (q = 0.793) (Table 8). (PRE:Precipitation; PET:Potential evapotranspiration).

Our results show that precipitation, soil type, potential evapotranspiration, altitude, and their interactions are the main reasons for the differences in the spatial distribution of WCA in the three basins of the study area.

When the precipitation exceeds the interception, depression detention, and infiltration of the underlying surface, surface runoff will occur (Tang et al. 2019). Regardless of whether it is a runoff generation under saturated conditions or a runoff resulting from excess rain, the surface runoff increases with the increase in precipitation, which further affects the WCF (Jiang et al. 2018). Evapotranspiration is the main form of surface water consumption affected by multiple factors such as atmospheric temperature, water vapor pressure, wind speed, solar radiation, and underlying surface conditions (Maček et al. 2018; Liu et al. 2019). Soil evaporation and vegetation transpiration are the primary forms of basin evapotranspiration, which together consume the water stored in the basin, and the amount of evapotranspiration directly affects the WCA of the basin (Li et al. 2016). Affected by various factors such as atmospheric circulation and topography, precipitation and evapotranspiration in mountainous areas have obvious spatial stratified heterogeneity, which affects the spatial distribution characteristics of WCA (Zhu et al. 2012).

The soil layer is the main body of terrestrial ecosystems for water conservation (Bormann 2012). Previous studies have shown that the smaller the soil bulk density, the relatively loose and porous the soil, the larger the capillary porosity, the faster the water seepage, making the soil layer store more water, thereby improving the soil's water conservation capacity (Bauer et al. 2014). As an inhomogeneous and changing space-time continuum, soil resources are affected by natural factors such as soil parent materials, topography, precipitation, and human activities, and have a high spatial stratified heterogeneity (Zhang et al. 2012; Reza et al. 2015). Furthermore, different soil types have different physical and chemical properties such as texture, organic matter, and soil bulk density. Therefore, the spatial stratified heterogeneity of soil significantly affects the spatial distribution of WCA.

As a natural geographic change, the altitude gradient directly affects the vertical distribution of precipitation and evapotranspiration, which affects the physical and chemical properties of soil, soil water and heat conditions, vegetation distribution, and vegetation density. It ultimately results in gradient changes in WCF (Li et al. 2012). The upstream of Xiong'an New Area is located in a mountainous area, and the topography in the area is very undulating. The elevation of different regions is quite different in this area, which leads to the difference in the spatial distribution of WCA.

The SWAT model constructed for the three catchments upstream of the Xiong'an New Area has high accuracy, and can reasonably reflect the hydrological cycle process at the annual and monthly scales for the three basins, and provides a tool for the accurate estimation of the basin WCA. At an annual scale, during the study period, the inter-annual changes in WCA in the ZJG, ZTM, and FP basins were large, and the average multi-year WCA was 47.59 mm, 44.59 mm, and 4.55 mm, respectively. The annual WCA of the ZTM basin was always positive. The WCA of the ZJG basin showed negative values in 2009 and 2014, whilst the FP basin showed negative values in 2009, 2013 and 2014. Spatially, the WCA in parts of the upstream of the ZJG basin, parts of the middle and upper reaches of the ZTM basin, and the FP basin has undergone a positive and negative conversion between years. The areas where the WCA did not change between positive and negative values did have changes in the high and low values of the WCA. On a monthly scale, there are apparent monthly differences in WCA in the ZJG, ZTM, and FP basins. The WCA of the rainy season months (June-September) are all positive, and the WCA of the non-rainy months (May-October) are all negative. Spatially, each sub-basin of the ZJG, ZTM, and FP basins has a positive and negative value conversion of WCA between the rainy season and the non-rainy season. Precipitation, evapotranspiration, and their combined effects are the main factors leading to the temporal changes and spatial pattern of WCA in the study area.

During the construction and operation management of the Xiong'an New Area, climate factors should be fully considered, and the urban functional layout and construction standards of significant infrastructure should be scientifically designed to reduce the adverse effects of climate change under the background of global warming.

This research was funded by the National Natural Science Foundation of China (No. 41877170), National Key Research and Development Program of China (No. 2018YFC0406501-02) and Hebei Province Key Research and Development Program of China (No. 20324203D and No. 20536001D). This work is supported by CFERN & BEIJING TECHNO SOLUTIONS Award Funds for excellent academic achievements.

Data cannot be made publicly available; readers should contact the corresponding author for details.