Regional water cycle response to land use/cover change for a typical agricultural area, North China Plain

Land use/cover (LULC) change is one of the most important drivers of ecohydrological changes (Blöschl et al. 2007; Li et al. 2009). Human economic activities represented by deforestation, reforestation, land reclamation, and urbanization alter land cover status dramatically and affect the surface hydrological processes (Legesse et al. 2003; Mao & Cherkauer 2009; Wijesekara et al. 2012; Lejeune et al. 2014; Bai et al. 2018) and land–atmosphere interaction (Chu et al. 2010; Huang & Margulis 2010; Zhang et al. 2013). LULC change not only affects global environmental changes (Gedney et al. 2006; Bonan 2008), but also has a significant impact on the regional water cycle (Ruprecht & Stoneman 1993), and thus affect climate change indirectly (Zhang et al. 2013).

The natural circulation of water near the surface of the Earth includes several important parts: precipitation, evapotranspiration, infiltration, surface runoff (SURQ), and groundwater (Bari & Smettem 2004; Liu et al. 2008; Zucco et al. 2014; Chawla & Mujumdar 2015; Meng & Wang 2017). Human economic activities affect every part of the hydrological cycle by changing LULC patterns (Allen & Ingram 2002; Sandra & Sathian 2016). Thus, the effects of such developments on components of the hydrological cycle assume significant importance, especially for the rational utilization of water resources in agricultural areas.

Many studies have investigated the relationships between LULC patterns and hydrological processes (Wang et al. 2012; Liu et al. 2013; Zhang et al. 2013; Cheng et al. 2017). The results show that hydrological processes are affected not only by LULC types (Jian et al. 2015; Duan et al. 2016) but also by spatiotemporal heterogeneity (Chu et al. 2010; Liu et al. 2013). However, debate on the effects of LULC changes still exists due to the complex interaction between human activities and the hydrological cycle (Lørup et al. 1998; Alkama et al. 2013; Liu et al. 2016). In addition, many studies have estimated the hydrological response to LULC changes at the basin scale (Zhang et al. 2014, 2016a, 2016b; Cheng et al. 2017). Less attention has been paid to the spatiotemporal heterogeneity of the LULC change effects on hydrological components at the sub-basin scale.

Generally, there are two methodologies to assess the impact of LULC change on hydrology. One is based on paired catchment experiments. This approach is time-consuming and difficult to apply to large basins (Zuo et al. 2016). The other is constructing hydrological models to explore the response of the hydrological cycle to LULC changes under different scenarios (Bao et al. 2019). Brown et al. (2013) analyzed the impacts of vegetation changes on streamflow regimes based on this method and reported that impacts on low flows were greater than those on high and/or median flows. Many studies also focused on the effects of vegetation on hydrological components (Bosch & Hewlett 1982; Zhang et al. 2011). Among all the studies, the responses of different hydrological components to LULC changes and its spatiotemporal features have rarely been considered (Brown et al. 2005; Wang et al. 2012; Jian et al. 2015).

Hydrological modeling is a useful tool for evaluating the response of different hydrological components to LULC change (Yang et al. 2017; Lin et al. 2018; Zhi et al. 2018). The Soil and Water Assessment Tool (SWAT) model has been widely used for its ability to analyze the spatiotemporal heterogeneity of influence by considering topography, soil, land-use and other basin characteristics (Allen & Ingram 2002; Sandra & Sathian 2016; Meng et al. 2018). The accuracy of hydrological simulation will be greatly limited for basins lacking meteorological data. The China Meteorological Assimilation Driving Datasets (CMADS) (version 1.0) developed by the China Institute of Water Resources and Hydropower Research (IWHR) (Meng et al. 2018) can solve this problem with higher resolution and better accuracy. The datasets range from 2008 to 2016 and cover the entire East Asian region (Meng & Wang 2017). Some studies considered that CMADS + SWAT has better results for runoff simulation (Liu et al. 2018).

In this study, a novel framework was proposed to quantify the effects of LULC change on the hydrological components and analyze the spatiotemporal heterogeneity of the influences. First, the spatial and temporal differences in LULC change were analyzed in detail. Second, CMADS + SWAT was applied for hydrological simulation because the only three meteorological stations had a poor representation of the meteorological conditions in the study area. The SWAT model was calibrated and validated using parameters associated with different land cover maps. Third, different LULC scenarios were used to simulate the several phases of the hydrological cycle based on the SWAT model in the study area. The influences of LULC change on key components of the hydrological process were analyzed at both basin and sub-basin scales.

Xiaoqing River is one of the major rivers in the central area of Shandong Province, China. It originates from Jinan City and flows through Licheng, Zhangqiu, Zouping, Guangrao, Shouguang, and then into Laizhou Bay, Bohai Sea. The areas it flows through are all important agricultural areas in China. The Xiaoqing River Basin (XQB) has an area of 10,366 km² (116°30′–118°45′E and 36°10′∼37°30′N), with the elevation ranging from −12 m to 1,046 m (Figure 1). Flat topography characterizes most of the basin except for the southern mountains. The study region is dominated by the temperate monsoon climate, with an average annual temperature and average annual precipitation of 13.4 °C and 646.7 mm. About 50–70% of the annual precipitation is concentrated in June to September, and rainstorms generally occur in July and August. The total population in this region exceeds 27.4 million in 2017, with Jinan, the capital of Shandong Province, accounts for 34.10% of the total (Li et al. 2020). The main land-use types in this basin are forest (FRST), urban area (URBN), and cropland (AGRL), followed by water body (WATR), saline land (WETL), and barren land (BARR). Among them, the total area of forest, urban and cropland exceeds 85% (Li et al. 2020) of the total. The main soil types in the XQB are cinnamon soil and moist soil.

The basic data used in this study are listed in Table 1 and described as follows:

1. The first version of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) (grid cell: 10 m × 10 m) was used for sub-basin division, river system extraction, and slope reclassification.

2. Land-use data for 2010 and 2016 were interpreted from two cloud-free Landsat 8 OLI images with a resolution of 30 m and support vector machine algorithm. The LULC types include WETL, WATR, URBN, BARR, FRST, and AGRL (Figure 2).

3. The soil input data were obtained from the 1:1 million soil dataset created by the Cold and Arid Regions Sciences Data Center at Lanzhou. The spatial resolution of this soil database is 1 km × 1 km. The soil data were resampled before spatial analysis.

4. The meteorological data were taken from the CMADS version 1.0. This dataset include precipitation, temperature, relative humidity, solar radiation, wind speed, location, and the elevation of each site. Temperature, relative humidity, and wind speed were generated from 2,421 national automatic stations and 39,439 regional automatic stations. Precipitation data were generated from the integration of multiple satellite data and precipitation from ground automatic stations. The production of radiation data was based on the Discrete Ordinates Radiative Transfer (DISORT) radiative transfer model and the acquisition of products from the FY2E satellite primary product for inversion of solar shortwave radiation. Twenty CMADS meteorological stations were used in this study (Figure 1).

5. The hydrological data were provided by the Nanjing Hydraulic Research Institute and comprised measured monthly data from 2008 to 2016 at the Shicun Hydrological Station.

Based on the land-use data generated from Landsat 8 OLI images, the LULC transfer matrix was calculated between 2010 and 2016. To analyze the spatial difference of LULC change, the conversions of land-use types were calculated at the sub-basin scale.

The SWAT model was developed by the U.S. Department of Agriculture (USDA) – Agricultural Research Services (ARS). It is a semi-distributed continuous hydrologic model (Arnold et al. 2012). The entire basin simulation domain is divided into sub-basins that are further subdivided into uniform hydrological response units (HRUs) with homogeneous soil, LULC, and slope characteristics. In this study, the XQB was divided into 29 sub-basins and 2,454 HRUs. The model parameters varies among different HRUs due to different land surface conditions. Theses parameters were determined through model calibration and validation. The simulation in this study was divided into three periods: warm-up period (2008–2009), calibration period (2010–2014), and validation period (2015–2016). Land-use data for 2010 were used for warm-up period and calibration period, while land-use data for 2016 were used for validation period.

Model calibration and parameter sensitivity analysis were conducted based on the SWAT Calibration and Uncertainty Program (SWAT-CUP). SWAT-CUP is a public domain program developed by Abbaspour (2007). The Sequential Uncertainty Fitting version 2 (SUFI-2) algorithm in SWAT-CUP was selected as the calibration algorithm due to the ideal performance in large basins (Rostamian et al. 2008). It provides a comprehensive optimization and uncertainty analysis through the global search method (Abbaspour 2014).

The SURQ was estimated by the modified Soil Conservation Service (SCS) Curve Number method (Neitsch et al. 2011). In this method, the amount of runoff was estimated based on local land-use types, soil types, and the antecedent soil moisture conditions. The Penman–Monteith method was used to estimate potential evapotranspiration effects. The physically meaningful absolute minimum and maximum ranges and fitted value for the key parameters involved in these methods were established (Table 2). Some studies focus on the analysis of parameter selection and sensitivity of the SWAT model (Cibin et al. 2010; Guse et al. 2013).

URBN, FRST, and AGRL were the primary land-use types in the XQB. The area of these three types accounted for 88 and 86% of the total area in 2010 and 2016, respectively (Table 3). The area of URBN and FRST showed upward trends totally. The area of URBN increased about 25.57% and its proportion of the total basin area has increased from 26% in 2010 to 32% in 2016. Similarly, the area of FRST increased about 10.56% and its proportion of the total basin area has increased from 30 to 34%. However, the area of AGRL decreased about 36.76% and its proportion of the total basin area has decreased significantly from 32 to 20%. The areas of WETL, WATR, and BARR remained relatively fixed in the study period.

To further analyze the change trajectory of main land-use types, the transfer matrices of URBN, FRST, and AGRL were calculated at the basin scale (Table 4). Overall, changes mainly occurred among these three types. About 23.57% of URBN converted to FRST and AGRL, and more than 19% of FRST converted to URBN and AGRL. Besides, nearly 50% of AGRL converted to URBN and FRST. These mainly be due to the relocation of towns in some southern mountainous areas and the implementation of ‘Grain and Green’ (Figure 2).

To identify the spatial and temporal differences of land-use change, the transferred matrix was also calculated at the sub-basin scale (Figure 3). The net conversion rate was calculated by the difference between conversion rates of two land-use types. The results showed that 6.78% of FRST and 10.27% of AGRL converted to URBN, respectively (Table 4). The sub-basins where AGRL/FRST converted to URBN were mainly distributed in the upper reaches (sub-basins 17, 20, 24, 27, and 28) and the lower reaches (sub-basins 1–5) of the whole basin. This is mainly due to the urbanization and population increase around Jinan City in the upper reaches and coastal areas recently (Li et al. 2020). In addition, 5.89% of AGRL converted to FRST in the XQB (Table 4). The conversion of AGRL to FRST occurred in most sub-basins (Figure 3). The results indicated that evident LULC changes may be due to the implementation of the ‘Grain for Green’ in the XQB.

The simulation result of monthly runoff from 2010 to 2016 in the XQB is shown in Figure 4. The R² and NSE for the calibration and validation periods were 0.81, 0.77 and 0.69, 0.60, respectively. Although the simulated result slightly underestimated streamflow in the dry period, the two statistical indexes (R² and NSE) indicated that the modeling accuracy was acceptable (Table 5).

Potential reasons for the underestimation can be categorized as measurement quality and model behavior. Measurement quality refers to missing precipitation data or streamflow measuring error. The study area is distributed in an agricultural plain region. Smaller slopes and larger amounts of river networks result in highly variable precipitation and difficulties in spatially estimating the precipitation (Jaepil et al. 2013). Additionally, considering the reduced accuracy caused by fewer meteorological stations in the XQB, the CMADS meteorological data were used for simulation (Figure 1). It has been reported that the accuracy of streamflow prediction mainly depends on the precipitation gage numbers and corresponding locations (Cao et al. 2006; Mul et al. 2009). Thus, the insufficient meteorological records and the distances among the meteorological stations result in limited simulation accuracy.

Model behavior refers to model limitations such as the inadequate representation of the physical mechanism of hydrological processes. For example, the SWAT uses total daily precipitation without considering rainfall intensity. Thus, runoff for some precipitation events (Qiu et al. 2012) can be estimated. Another example is the use of runoff curve numbers to simulate the SURQ behavior. This approach does not account for saturation excess runoff or contributions from variable source areas (Garen & Moore 2005; Easton et al. 2008).

Monthly runoff from 2010 to 2016 was simulated based on two LULC scenarios (2010LULC and 2016LULC) using the adjusted model. Three hydrological component indicators were analyzed: SURQ, ET, and WYLD. The total average values of these hydrological components were calculated (Figure 5(a)). SURQ increased by 49% and WYLD decreased by 47%, while ET remained relatively stable. The multi-year monthly average values of the hydrological components were also calculated (Figure 5(b)–5(d)). The hydrological components showed obvious seasonal variations, with high change values during July and August. SURQ increased 2.47–2.72 × 10⁴ mm in these two months under two LULC scenarios (Figure 5(b)) and WYLD increased 1.61–2.03 × 10⁴mm (Figure 5(d)). The multi-year monthly average values of SURQ and WYLD under 2016LULC were 1.0–1.6 and 0.9–2.0 times of that under 2010LULC, respectively. However, there was little difference between the results of simulated ET under two scenarios (Figure 5(c)).

To further analyze the influences of different land-use types on the hydrological components, the proportions of hydrological components produced by different land-use types to the total amount of hydrological components were calculated (Table 6). URBN was the main type that generated SURQ under two scenarios. About half of the SURQ in the basin was generated from URBN. The expansion of the impervious surfaces caused by urbanization prevents water infiltration and brings more SURQ. From 2010 to 2016, the area of FRST increased from 3,142 to 3,473 km² (Table 3). The proportion of SURQ generated from FRST in the total amount of SURQ generated by all land-use types under two scenarios had decreased by 2.91%. The expansion of forest area has improved the capacity to conserve water in the study area.

Moreover, WATR, AGRL, and FRST were the main land-use types that generate ET. From 2010 to 2016, the proportions of ET generated on AGRL and FRST in the total amount of ET generated by all land-use types had increased by 0.89 and 0.55%, respectively. Since the beginning of the 21st century, the policy of ‘Grain for Green’ has been gradually implemented, especially in Northern China (Schilling et al. 2018). Forests can store more water than cropland, and decrease streamflow and increase ET (Schilling et al. 2018) simultaneously. Meanwhile, about 5.89% of AGRL converted into FRST (Table 4). Therefore, under the same precipitation, more ET could be generated as the percentage of FRST increased from 2010 to 2016.

URBN, FRST, and AGRL were the main land-use types that generated WYLD. From 2010 to 2016, the WYLD generated on URBN in the total amount of WYLD generated by all land-use types had increased by 9.55%, with other land-use types showed decreasing trends. The influence mechanism of land-use on WYLD is complicated. In the SWAT model, WYLD is also related to SURQ, lateral flow, groundwater, and transmission losses in the shallow aquifer (Arnold et al. 1993). AGRL could hold the water in plants and soil. Meanwhile, crops need water to grow. FRST could transport water into the soil and produce lower WYLD (Im et al. 2003; Yang et al. 2013). In this study, the WYLD generated on FRST in the total amount of WYLD generated by all land-use types had decreased by 2.34% from 2010 to 2016, with the percentage of FRST area has risen by 4% (Table 3) in the same period. This was consistent with most previous studies (Bi et al. 2009), which suggested that forest may be an important factor for the WYLD.

The spatial distribution of changes of ET, SURQ, and WYLD at the sub-basin scale is shown in Figure 6. SURQ increased more than 15% in some central and downstream parts of the whole basin, while minor increases occurred in some upstream sub-basins (e.g. sub-basins 17 and 20). ET increased obviously in the southern part and some upstream areas. The areas that were closed to the outlet also showed a 0.02–0.05% increase (e.g. sub-basins 2 and 3). As for WYLD, only some upstream parts (e.g. sub-basins 15 and 17) and the areas near the estuary (e.g. sub-basins 1 and 2) showed a significant increase, while other areas showed a downward trend, especially the southern parts of the whole basin.

Quantitative statistical analysis of the impacts of LULC changes on these hydrological components at the sub-basin scale is shown in Figure 7. The results indicated that the hydrological components varied among different LULC types in the whole basin.

URBN has the largest contribution rate to SURQ in the downstream part near the estuary and upstream part of the basin. Taking the sub-basin 20 close to Jinan City as an example, the SURQ generated on URBN has increased by more than 11%. A similar trend also occurred in sub-basin 4 (Figure 7). Numerous studies showed that urbanization alters a basin's response to precipitation events, leading to increased volumes of SURQ (McColl & Aggett 2007; Oudin et al. 2018; Li et al. 2020). In this study, the impervious areas expanded obviously due to the urbanization in the coastal areas and upstream areas. It leads to the weakening of the infiltration capacity of the ground and increasing the percentage of impervious surfaces (Cuo et al. 2008; Wang et al. 2012).

In contrast, ET showed a slightly decreasing trend in most areas of the XQB (Figure 6(c)). For example, ET generated on AGRL in sub-basin 18 decreased by 5.61% (Figure 7), where AGRL was the main land-use type. A similar trend also appeared in regions such as sub-basins 12 and 13. However, for FRST-dominated areas, the amount of ET increased as the area of FRST expanded. For example, the ET generated on FRST in sub-basin 29 increased by 6.53% and more than 15% of AGRL converted to FRST in the same sub-basin (Figure 3). Other studies reveal that ET rates from FRST are higher than other LULC types, generally due to their canopy interception and transpiration rates (Guse et al. 2013; Locatelli et al. 2017). In this study, the increasing trend of FRST in some sub-basins may be due to the implementation of the ‘Grain for Green’ policy that was mentioned above.

Numerous studies show that the combined effect of urbanization and forest reduction leads to an increase in WYLD, which may be due to an increase in direct runoff from impervious surfaces (Zhang et al. 2016a; Oudin et al. 2018). In the current study, the sum of WYLD generated on URBN and FRST increased about 10.43% in sub-basin 4, while that on AGRL decreased about 2.61% (Figure 7). From the perspective of the LULC transfer matrix (Figure 3), the combined increase of URBN and AGRL exceeded 20%, while the FRST decreased about 10% in the same region. Similar situations also can be seen in sub-basins 8 and 20 where urbanization and deforestation were obvious (Figure 3). These findings were consistent with previous studies. Nevertheless, the process of WYLD was complicated, especially when other hydrological components changed, such as groundwater and lateral flow (Li et al. 2018). Previous studies showed that the change of WYLD is also related to such as climate change, regional differences, and soil conditions (Im et al. 2003; Li et al. 2018). Although the results of this study can reflect the general tendency of WYLD and reveal the relationship between WYLD and LULC to some degree, it is still necessary to further investigate in the future.

This study quantified the impacts of spatiotemporal LULC change on the water cycle at both basin and sub-basin scales in the XQB, a typical agricultural area located in the North China Plain. The SWAT model was used for monthly runoff simulation under two LULC scenarios. Both R² and NSE reached acceptable values in the calibration and validation periods. URBN, FRST, and AGRL were the dominant LULC types in the XQB. The land-use patterns changed significantly from 2010 to 2016. The urban and forest areas increased by 25.57 and 10.56%, with cropland area decreased by 36.76% from 2010 to 2016. From the perspective of net conversion rates, about 6.78% of FRST and 10.27% of AGRL converted to URBN, respectively. Besides, about 5.89% of AGRL converted to FRST in the XQB. These trends were closely associated with the development of urbanization and the implementation of the ‘Grain for Green’.

The impacts of LULC change on hydrological components were also assessed quantitatively. Our analysis showed that SURQ has increased by 49% due to the expanded impervious surface area (25.57%) caused by the development of urbanization, especially in some upper and lower sub-basins. In other words, runoff increases as permanent imperviousness increases. There was little change of ET (−0.03%) under two LULC scenarios at the basin scale. However, ET increased in some sub-basins where the conversion of AGRL to FRST was obvious, which was the result of the implementation of ‘Grain for Green’. For instance, ET increased by 6.53% in sub-basin 29, with 15% of the cropland area converted into forest. WYLD has increased by 47% from 2010 to 2016 at the basin scale generally. We roughly revealed that the increase of WYLD in some sub-basins was related to urbanization and forest reduction. These differences can explain the relationships between hydrological characteristics and LULC change at both basin and sub-basin scales to some extent.

In this study, the impacts of LULC change on the model performances were analyzed based on the calibrated results and uncertainty analysis of the model outputs. The significance of using dynamic land-use input instead of static land-use input when simulating hydrological components in a basin was revealed. This study can provide useful information for LULC planning and soil–water conservation in the agricultural area, North China Plain.

This study is funded by the National Key Research and Development Program of China (2018YFC1407601) and appreciates valuable comments and suggestions from editors and anonymous reviewers.

All relevant data are included in the paper or its Supplementary Information.