Meteorological data, soil data, and DEM are the input data for the model simulation. Hydrologic data are used to verify the simulated effect. The rationality and accuracy of the input data are all very important for the simulation performance of the model.

3.1.1 Meteorological data

The meteorological database comprises daily measurements collected at Anshun weather station by the China Meteorological Data Service Center (http://cdc.cma.gov.cn/home.do). The database includes daily average temperature, precipitation, average wind speed, relative humidity, and sunshine duration. In addition, because this study focuses on runoff simulations, which are more prone to the influence of precipitation, four precipitation stations in the Sancha River Basin were selected for precipitation data. These stations are located at Bide, Machang, Santang, and Qibo. Measurements at the precipitation stations were collected from the Institute of Geographic Sciences and Natural Resources Research, Chinese Academy of Sciences (CAS): Hydrological Yearbook of the People’s Republic of China - the Wujiang River Region of Yangtze River Basin. Basic information from weather and precipitation stations is shown in Table 1.

3.1.2 Hydrological data

According to the Yangtze River Basin record in the Hydrological Yearbook of the People’s Republic of China (Ministry of Water Resources of the People’s Republic of China), there are two hydrological stations in the Sancha River Basin: the Yangchang station and the Longchangqiao station, which are located in the upstream and midstream regions, respectively. Runoff data from hydrological stations in the discharge areas of the basin are required for model calibration, so Hongjiadu station (at the Liuchong River mouth) and Yachi station (in the mainstream of the Wujiang River, where the Liuchong River and Sancha River converge) were selected. Daily and monthly average runoff data from the Yangchang, Longchangqiao, Hongjiadu, and Yachi stations were used to calibrate and validate the SWAT model. Information about these stations is listed in Table 1.

3.1.3 Remote sensing data

Remote sensing data include digital elevation model (DEM) data, land use data, and soil data, which include soil types and associated physical and chemical properties. The Normalized Difference Vegetation Index (NDVI) was adopted for the spatial statistical analysis of runoff model simulation results. DEM data at 1 km resolution were obtained from the HYDRO1K database at the U.S. Geological Survey (USGS) EROS Center (http://eros.usgs. ov/#/Find_Data/Products_and_Data_Available/HYDRO1K). The 1:100,000 scale land use vector data of 2010 obtained from the Resources and Environmental Science Data Center, CAS (Liu et al., 2014) were applied. Soil data were collected from the Harmonized World Database version 1.1 (HWSD), established jointly by the Food and Agriculture Organization of the United Nations (FAO) and the International Institute for Applied System Analysis (IIASA) (Fischer et al., 2008). The main soil properties include soil type according to the FAO90 soil classification scheme, reference depth, and physical (such as gravel volume percentage, sand/silt/clay content, and effective moisture content) and chemical properties (namely organic carbon content and cation exchange capacity). NDVI data and the Moderate Resolution Imaging Spectroradiometer (MODIS) product (http://ladsweb.ascom. asa.gov/) from 2005-2010 were used, with temporal and spatial resolutions of 16 d and 500 m, respectively. Vegetation index data were mainly used to study the statistical spatial relationship between land cover and runoff generation in the river basin. NDVI in the Sancha River Basin is higher in the west and lower in the east, showing a banded spatial distribution (Figure 2).

This study used the SWAT model, which was calibrated and validated using measured runoff data from hydrological stations. It was then utilized to simulate runoff generation in a typical karst river basin, including total runoff, surface runoff, and groundwater runoff. The water cycle in the model was based on various hydrological variables:

$S{{W}_{t}}=S{{W}_{0}}+\sum\limits_{i=1}^{t}{({{R}_{day}}-{{Q}_{surf}}-{{E}_{a}}-{{W}_{seep}}-{{Q}_{gw}})}$ (1)

where SW_t denotes the final soil moisture content (mm); SW₀ denotes the initial soil moisture content on the i th day (mm); t represents the time (d); R_day represents the precipitation on the ith day (mm); Q_surf is the surface runoff on the ith day (mm); E_a is the evapotranspiration on the i the day (mm); W_seep denotes the water seepage from the soil profile to the vadose zone (mm); and Q_gw denotes the baseflow returned on the i th day (mm).

Compared to other hydrological models (distributed and conceptual hydrologic models), the SWAT model is more favorable because it is based on hydrological processes and is more efficient due to easily obtained parameters. For example, the MIKE SHE hydrological model comprehensively considers hydrological processes, boundary conditions, and the spatial heterogeneity of collective characteristics in river basins; however, it is difficult to collect data on underlying surfaces and hydrogeological conditions (McMichael et al., 2006). The TOPMODEL model has a simple structure and demands fewer parameters, but it does not thoroughly consider the spatial heterogeneity of hydrological elements, and the coupled relationships between hydrological units (Beven et al., 1984). In the SWAT model, runoff simulation is performed using different modules (surface runoff, groundwater, and river confluence). For surface runoff, this study employs the relatively mature SCS curve number method (empirical model), which is a function of soil permeability, land use, and antecedent soil moisture conditions. The soil permeability is predominantly calculated according to the physical and chemical properties of different soil types in the karst area, which can reasonably reflect actual soil permeability in the region. Different land use types will affect runoff by changing surface evaporation, soil moisture conditions, and interception of land cover. As an input data point, the land use types in the river basin are eventually reflected by the SCS curve number. Studies have indicated that surface runoff will occur after the precipitation depth reaches 40 mm in karst areas (Peng and Wang, 2012). Therefore, antecedent soil moisture conditions are also a key factor influencing surface runoff in typical karst river basins. Groundwater runoff was calculated using water mass balance equations of shallow and deep aquifers. The volume was calculated based on different movement processes, and each module was controlled by the corresponding parameters. Groundwater runoff simulation results in karst areas are generally affected by factors such as the baseflow alpha factor (ALPHA_BF, a direct indicator of groundwater runoff in response to recharge variations) and threshold water depth in the shallow aquifer for flow (GWQMN).

3.2.1 Data preparation and model construction

(1) Construction of spatial database

Coordinate systems of different spatial data were consolidated (Krasovsky_1940_Albers) through ArcGIS using the .img data format. According to the land use data classification in the SWAT model, land use raster data were re-classified and converted to codes defined by the model. The 1:1,000,000 soil data for the Sancha River Basin underwent projection conversion, and special situations in the attribute tables were determined according to previous research experience.

(2) Preparation of attribute database

Soil database. The major parameters required by “USERSOIL” in the SWAT database are soil physical and chemical properties. Physical properties include mechanical composition, saturated hydraulic conductivity, and available water capacity, and chemical properties include organic carbon content, electrical conductivity, and others. Some parameters were obtained from the HWSD soil database, for example, soil layer thickness and organic carbon content, and saturated hydraulic conductivity and available water content were calculated using software (SPAW) or empirical equations.

Meteorological database. For each weather station, the .DBF documents of daily precipitation, daily maximum and minimum temperatures, solar radiation, average wind speed, and relative humidity were established. Meanwhile, a geographic information document (.DBF) of weather and precipitation stations was created to link the location map of the weather stations with the meteorological database.

Land use and soil index table. The index table links the SWAT attribute data with the spatial distribution maps of land use and soil type.

3.2.2 Sensitivity analysis and parameter calibration of SWAT model

This study utilizes measured runoff data from the Yangchang, Longchangqiao, Yachi, and Hongjiadu hydrological stations for parameter sensitivity analyses and model calibration and validation. Parameter sensitivity analyses were performed using the LH-OAT method embedded in the model and measured runoff data from hydrological stations. The results are shown in Table 2, which demonstrates that the most influential parameters for the runoff simulation were: the initial curve number (II) value (CN₂), which controls surface hydrological processes under soil moisture condition II (moderately moist soils); available water capacity (SOL_AWC); the soil evaporation compensation factor (ESCO); maximum canopy storage (CANMX); baseflow alpha factor (ALPHA_BF), which controls subsurface hydrological processes; the groundwater delay time (GW_DELAY); threshold water depth in the shallow aquifer for flow(GWQMN); and channel effective hydraulic conductivity (CH_K₂), which controls the main stream confluence. Based on the parameter sensitivity results and actual hydrological processes, the adjustable parameters were CN₂, SOL_AWC, ALPHA_BF, GWQMN, and CH_K₂. These parameters were adjusted several times until the simulated values were close to the measured values. The final values for all parameters are listed in Table 2. For groundwater runoff, previous studies on subsurface hydrological processes in karst areas (recharge, runoff, and discharge) (Amatya et al., 2011; Chen et al., 2014) were examined. Parameters were adjusted and validated several times to ensure accurate simulation results. Eventually, ALPHA_BF and GWQMN were set to 0.048-0.5 and 1000, respectively.

The Nash-Sutcliffe efficiency, E_ns, coefficient of determination, r², and hydrograph based on measured monthly runoff data for 2008-2010 were applied to evaluate the efficiency of the model and validate the calibration results. E_ns values closer to 1 indicate simulation values that are closer to the measured values. E_ns was calculated as follows:

${{E}_{ns}}=1-\frac{\sum\limits_{i=1}^{n}{{{({{O}_{i}}-{{P}_{i}})}^{2}}}}{\sum\limits_{i=1}^{n}{{{({{O}_{i}}-\overline{O})}^{2}}}}$ (2)

where O_i, O, P_i, and P are measured values, average measured values, simulated values, and average simulated values, respectively. It is generally believed that E_ns within 0.5-0.65 indicates acceptable simulation results. E_ns values of 0.65-0.75 and above 0.75 signify relatively good and excellent simulation results, respectively (Popov, 1979).

The results suggest that simulated monthly average precipitation at Longchangqiao station matches the measured values relatively well during the calibration period, with an E_ns of 0.82 and an R² of 0.92 (Figure 3 and Table 3). During the validation period, E_ns and R² are 0.90 and 0.93, respectively. The simulation agrees well with the measurement at Yangchang station during both calibration and validation periods, with corresponding E_ns values of 0.70 and 0.73. The simulation results are thus relatively good. In addition, from the hydrograph, the simulated and measured values at Yangchang and Longchangqiao stations are relatively coherent. Furthermore, simulation results at the basin outlet of both stations are not as good as those of Longchangqiao and Yangchang stations. The measured runoff in the discharge areas is estimated based on the difference between Yachi and Hongjiadu stations, which may contribute to the partial discrepancy.

This study used the NDVI vegetation index to examine the spatial heterogeneity of runoff generation, and the land use landscape fragmentation index to explain the spatial relationship between runoff service spatial heterogeneity and runoff volumes. The degree of landscape fragmentation was measured using the effective mesh size (Jaeger, 2000), which refers to the average value of the continuous area of each land use type in a landscape. A small value denotes a less fragmented landscape. The landscape fragmentation index can spatially reflect the intensity of land use type variations among different hydrological response units (HRUs). In more fragmented areas, there are remarkable differences in land use types between neighboring HRUs. The index integrates ecological processes, landscape compositions, and spatial patterns, and can therefore comprehensively and objectively express the degree of fragmentation (Gao and Cai, 2010). The equation is as follows:

${{m}_{eff}}(j)={{A}_{j}}\cdot \sum\limits_{i=1}^{n}{{{(\frac{{{A}_{ij}}}{{{A}_{j}}})}^{2}}}=\frac{1}{{{A}_{j}}}\sum\limits_{i=1}^{n}{{{A}_{ij}}^{2}}$ (3)

where m(j), n, A_ij, and A_j denote the effective mesh size of landscape j, the number of complete patches in landscape j, the area of patch i in landscape j, and the area of landscape j, respectively. If the minimum of the index equals the size of the grid, it indicates that land use types are different for all neighboring grids. If the maximum equals the landscape area, it suggests that the landscape is of a homogeneous type. This study used the ‘moving window’ analysis (research amplitude) in the landscape pattern analyzing software FRAGSTATS (McGarigal and Marks, 1995) to compute the effective mesh size.

Based on the 30-m land use data (2010) in the Sancha River Basin (Figure 2b), it demonstrates that areas with high and low landscape fragmentation indices are scattered all over the region. The total area of low-value (highly fragmented) regions is larger than that of high-value regions. High-value areas are mainly located along the upstream region with higher elevations and less anthropogenic interference. Low-value areas are concentrated in the mid- and downstream regions.

Geographically Weighted Regression (GWR), proposed by Brunsdon et al. (1996; 1998), is a simple and practical local spatial analysis technique. GWR depicts the internal spatial relationships in the study area. In this study, it was used to reveal the spatial relationships of runoff volume with landscape fragmentation indices and the NDVI vegetation index under linear non-stationary conditions. The GWR model is an extension of the traditional regression model. It uses the spatial locations of data as one of the parameters, and evaluates spatial variations in independent variable-dependent variable relationships through using parameters.

The basic form of the GWR model is:

${{y}_{i}}={{\beta }_{0}}({{\mu }_{i}},{{\nu }_{i}})+\sum\limits_{k=1}^{p}{{{\beta }_{k}}({{\mu }_{i}},{{\nu }_{i}})}{{x}_{ik}}+{{\varepsilon }_{i}}$ (4)

where (μ_i, v_i) and k represent the coordinates of the ith sampling point and the number of independent variables, respectively; y_i, x_ik, and ε_i are the independent variables (total, surface, and groundwater runoff volumes), dependent variables (landscape fragmentation index and NDVI), and random errors at regression point i, respectively. β₀(μ_i, v_i) and β_k(μ_i, v_i) denote the intercept and slope of the GWR model at regression point i, respectively. The parameters can be estimated using the following formula:

$\beta ({{\mu }_{i}},{{\nu }_{i}})={{({{X}^{T}}W({{\mu }_{i}},{{\nu }_{i}})X)}^{-1}}{{X}^{T}}W({{\mu }_{i}},{{\nu }_{i}})Y$ (5)

where β(μ_i, v_i) and W(μ_i, v_i) are, respectively, the unbiased estimator of the regression coefficient and the spatial weighting matrix, which is the core of the GWR model; the weight increases when the observation point is closer to the defined point. X and Y are the matrices for independent and dependent variables, respectively.

The Gaussian Spatial Weight function is universal, and its expression formula is as follows:

${{\omega }_{ij}}=\exp \left( -\frac{d_{ij}^{2}}{{{b}^{2}}} \right)$ (6)

where ω_ij, d_ij, and b denote the weight of observation point j, the Euclidean distance between regression point i and observation point j, and the non-negative attenuation parameter describing the weight-distance function relationship (also called bandwidth), respectively. When the distance between observation points is greater than b, the weight rapidly approaches 0.

The spatial patterns of runoff generation in the Sancha River Basin, including total runoff volume (TOTAQ), surface runoff (SURQ), and groundwater runoff (GWQ), are obtained using the SWAT model. In addition, because evapotranspiration significantly affects water balance, this study also analyzes the simulation results of actual evapotranspiration (ET). The average total runoff volume over the study period is higher in the north and lower in the south of the Sancha River Basin (Figure 4a), and fluctuates within 510-1177 mm, with an average of 826.4 mm. The total runoff coefficient reaches 70%. Areas with higher total runoff volumes (900-1100 mm) are mostly located in the east and the north, and cover 22% of the basin. The total runoff volume in the south of the basin is smaller and mostly within 500-700 mm.

The overall surface runoff is relatively low in the river basin, with an average value of 276 mm, and shows remarkable spatial heterogeneity. The surface runoff coefficient is 23.9%. Surface runoff in 50% of the river basin is between 0 and 300 mm, and mainly located in the southern portion. Areas with more than 300 mm runoff are scattered in the east and the north. This is because the Sancha River Basin is situated in a typical karst area, and surface runoff patterns are highly heterogeneous. The unique surface-subsurface binary hydrological structures cause a considerable amount of surface water loss to the subsurface. Hence, there is relatively abundant subsurface runoff in the river basin (Figure 4c). Approximately 60%-70% of the region has an underground runoff volume of 500-700 mm. Underground runoff volume in upstream areas is typically over 500 mm.

The actual evapotranspiration amounts show no significant differences over the river basin. Most areas in upstream and downstream regions experience evapotranspiration of 100-300 mm (Figure 4d), covering 53% of the total area. Evapotranspiration is 300-500 mm in certain parts of central and southern regions. Compared to the spatial patterns of different runoff types, the distribution of actual evapotranspiration is relatively homogeneous in the Sancha River Basin. This is because evapotranspiration is mostly influenced by temperature, precipitation, wind speed, and vegetation type (Borba et al., 2012). The DEM data (Figure 1) show 70% of the area is within the range of 900 to 1500 m and no remarkable differences in temperature and wind speeds over the study area. Coniferous forests and shrubs are the dominant vegetation types (more than 80%), which results in relatively low spatial heterogeneity of evapotranspiration.

Elevation in the Sancha River Basin ranges between 911 and 2330 m. The corresponding average values and coefficients of variation for the total runoff volume, surface runoff, underground runoff, and actual evapotranspiration at elevation gradients of 911-1000 m, 1000-1500 m, and 1500-2330 m are obtained. Figure 4 and Table 4 show that variations in the total runoff volume with elevation are small (803-882 mm), and fluctuations in the coefficient of variation are insignificant (0.09-0.21), while surface runoff variations are relatively large. At elevations of 911-1000 m, the surface runoff is 401 mm and the coefficient of variation is 0.52. As the elevation increases, the surface runoff decreases to 154 mm and the coefficient of variation increases to 0.85. The average value of underground runoff gradually increases at higher elevation. The actual evapotranspiration exhibits no difference at different elevation gradients. In general, elevation has significant effects on surface and underground runoff in the Sancha River Basin. The percentage of underground runoff increases with elevation, while that of surface runoff decreases.

The gradient analysis of slope (Table 4) indicates that the total runoff volume increases on steeper slopes. The surface runoff first increases and then decreases with increasing slope. It reaches a maximum (319 mm) on slopes from 5-10°. For a slope range of 0-10°, the surface runoff increases with slope. The force in the water layer along the slope surface increases, and the pressure perpendicular to the slope surface decreases. This leads to a smaller infiltration rate and greater surface runoff. As the slope becomes steeper, the vertical action force exerted by precipitation on the ground surface is smaller, and soil crust formation is slower, resulting in less surface runoff. This is similar to the law of runoff generation on karst slopes studied by Hu et al. (2012), whereby groundwater runoff essentially increases on steeper slopes. The highest runoff value is 629 mm. For steeper slopes, the vertical action force exerted by precipitation on the ground surface is smaller, and underground runoff increases at slower rates (Wu and Zhang, 2006). In different slope ranges, the average actual evapotranspiration fluctuates between 222 and 259 mm.

4.3.1 Statistical characteristics of runoff services under different land cover conditions

Vegetation cover influences runoff by reducing the kinetic energy of raindrops, intercepting precipitation, modifying soil crusts, and altering soil infiltration capacity (Wu and Zhang, 2006). The total runoff volumes are similar for the five land use types, and all fluctuations are within 718-850 mm. Total runoff volumes for grassland are slightly lower. There are relatively significant differences in surface and underground runoff among the different land use types. In particular, industrial and commercial land experiences the highest surface runoff (544 mm), followed by farmland. Commercial and industrial land is dominated by heavy human activities and impermeable layers in urban areas, which substantially weaken the capacity of soil to retain and store water. Consequently, most precipitation becomes surface runoff. Meanwhile, underground runoff shows the opposite trend. It accounts for more than 70% of the total runoff in orchard and forest land underground, and is the highest in forests (628 mm). Soils in forests have better physical structures, higher permeability, and the litterfall excels at retaining water (Shi and Li, 2001). Moreover, forests are often located in higher-altitude and steeper regions.

Vegetation cover conditions affect the hydrological processes and runoff generation in the river basin by changing surface evaporation, soil moisture conditions, and interception of land cover. The NDVI in the Sancha River Basin is divided into three classes: <0.5, 0.5-0.7, and >0.7, and the corresponding runoff for each NDVI gradient is analyzed. Figure 5 shows that the total runoff and underground runoff increase with higher NDVI, while the surface runoff increases at first but subsequently decreases (259-319-282 mm). This is because improved vegetation cover conditions reduce with limits resulting in water loss, which is due to surface-subsurface binary hydrological structures. As NDVI continues to rise, interception by the vegetation canopy increases and transpiration becomes stronger, and the soils can retain and store more water, which decreases surface runoff.

4.3.2 Spatial correlation between landscape fragmentation, vegetation cover, and runoff generation

Based on the GWR model, the landscape fragmentation index (m_eff), and NDVI, the local spatial relationships between land cover and different runoff volume variations are analyzed at the HRU scale (Figure 6). To ensure accurate characterization of macroscopic patterns and local details over the entire study area, R² and residual sum of squares of regression models are employed as criteria to determine the optimal bandwidth, which is 7 km. The Gaussian function is selected as the weight function. The total runoff volume and the landscape fragmentation index in the Sancha River Basin show mostly positive correlations. Areas with positive correlations occupy 96.5% of the entire region. Negative correlations mainly occur in limited areas in the southeast. Areas with strong positive spatial correlations are located in the south and northeast. The spatial regression relationship between underground runoff volume and the landscape fragmentation index is similar to that of the total runoff volume. 99.9% of the river basin shows positive correlations, and correlations gradually weaken from west to east. A relatively heterogeneous spatial relationship is noted between surface runoff and the landscape fragmentation index; they are mostly negatively correlated (74.2%) in western and central regions. Positive correlations are only observed in limited areas in the east and north. The actual evapotranspiration volume and the landscape fragmentation index demonstrate negative spatial correlations overall (96.9%).

The spatial non-stationary relationship of the GWR regression coefficient between NDVI and runoff shows remarkable differences in the south and west (Figure 7). In the upstream and midstream regions of the Sancha River Basin, the total runoff volume and NDVI are mostly positively correlated, which accounts for 61.4% of the entire basin. The correlations become negative in the north and the downstream region. The spatial relationship between surface runoff and NDVI is relatively complex, and the positively and negatively correlated areas are scattered throughout the basin. The positively correlated areas cover 67.3% of the entire region, and are mostly located in the western and central regions as well as the outlet of the basin. The underground runoff demonstrates negative correlations with NDVI for 56.4% of the basin. Figures 6 and 7 jointly suggest that areas with relatively better ecological background in the upstream region have higher elevation, less landscape fragmentation, greater NDVI, higher total runoff volume, more underground runoff, and less surface runoff. This indicates relatively effective ecological restoration. 52.5% of the basin shows negative correlations between actual evapotranspiration and NDVI, predominantly in the south and upstream regions.