Identifying driving factors of the runoff coefficient based on the geographic detector model in the upper reaches of Huaihe River Basin

2.1 Study area description

The study area is located in the upper reaches of Huaihe River Basin, including four major river basins: Shayin River, Hongru River, Huaihe River mainstream, and Shipi River. A total of 29 watersheds from the 4 major river basins with an area of between 50 and 2,000 km² were available for this investigation (Figure 1). The study area represents a transition zone between the northern subtropical region and the warm temperate zone. Rainfall and temperature progressively decrease from south to north. The mean annual precipitation ranges from 624 to 1,461 mm, with a mean value of 1,056 mm. The mean annual temperature ranges from 13.2 to 16.1°C, with a mean value of 15.1°C. The mean annual RC ranges from 0.17 to 0.68, with a mean value of 0.37. Of the four river systems, Shipi River has the highest values of mean annual precipitation, mean annual temperature, and mean annual RC. The characteristics of the 29 watersheds are given in Table 1. Mountainous areas with forest and grass are distributed mainly in the south and west. Farmland can be found in plain areas.

Figure 1

Hydrological stations of the study area.

2.2 Data sources

In this study, the runoff refers to the surface river. The mean annual RC for the 29 watersheds from 1980 to 2010 were obtained from hydrological yearbooks by the Ministry of Water Resources of China, and designated as the actual RC.

A total of 13 factors as the driving factors from 5 aspects of climate, topography, soil, vegetation, and land use type were representative and easy to quantify, and these data are readily available (Table 2). Precipitation, temperature, and ET are the most important climate factors. Precipitation and temperature datasets could be freely downloaded from the Resource and Environment Data Cloud Platform, Chinese Academy of Sciences (http://www.resdc.cn) with a 1,000 m resolution. Annual precipitation and temperature data were generated through spatial interpolation based on daily national meteorological observation data. The yearly 500 m MOD16A3 dataset from the Numerical Terradynamic Simulation Group at the University of Montana (http://www.ntsg.umt.edu) was applied to obtain actual evapotranspiration (AET) and potential evapotranspiration (PET) data. MOD16A3 dataset is based on the algorithm proposed by Mu et al. [11], which is based on the Penman–Monteith equation. The spatial distributions of these annual average precipitation, temperature, AET, and PET from 2001 to 2015 are shown in Figure 2a–e. The Shuttle Radar Topography Mission digital elevation model (SRTM DEM, version 4.1) with a 90 m resolution was used to extract elevation, slope, and topographic wetness index (TWI) (Figure 2e–g). The Harmonized World Soil Database was used to extract the sand, silt, and clay fractions of the topsoil [34] (Figure 2h–j), which was derived from the International Institute for Applied Systems Analysis (IIASA) (https://iiasa.ac.at/models-and-data/harmonized-world-soil-database). Saturated hydraulic conductivity (Figure 2k) was calculated using multiple pedotransfer functions [35]. The annual average normalized difference vegetation index (NDVI) was derived from SPOT VEGETATION NDVI time series datasets between 2001 and 2015 (Figure 2l). The land use map (Figure 2m) with 1 km resolution in 2010 included 6 categories: cropland, forest, grassland, water area, construction land, and unused land. The area percentages for each land use type of each watershed were computed from the land use map. The three main land use types (cropland, forest, and construction land) were selected.

Figure 2

Spatial distributions of 13 driving factors: precipitation (a), temperature (b), AET (c), PET (d), elevation (e), slope (f), TWI (g), sand (h), clay (i), silt (j), K _S (k), NDVI (l), and land use (m).

2.3 Water balance equation

For a given basin, the mean annual water balance equation is represented as:

where P, R, and AET are mean annual precipitation, runoff, and AET, respectively. ΔS is the mean annual change in water storage. Over a long period, the amount of water storage is zero. So, RC is represented as

The RC map can be estimated based on the precipitation map and the AET map from remote sensing data.

2.4 Geographical detector

Geographical detector proposed by Wang et al. [29] is a useful statistical method that is not only able to quantitatively detect the relationship between independent variables X and dependent variables Y but it can also be used to investigate the interaction between two independent variables to a dependent variable. Any assumptions of linearity and a specific form of interaction should not be considered. In this study, the factor detector and the interaction detector were adopted to explore the driving factors of the RC in the upper reaches of the Huaihe River Basin. The dependent variable Y was RC, while the independent variable X represented the 13 potential driving factors.

The purpose of the factor detector is to quantitatively detect the explanatory powers of different driving factors of the RC using the q-value. The formula is expressed as

where h = 1, 2,… L is the stratum of factor X (or the variables Y); N and N _h are the unit numbers of cells in the whole region and the layer h, respectively; and σ 2 and σ h 2 are the variances of variables Y of the whole region and the layer h, respectively. q ∈ [0,1]. A larger q value indicates that Y can be more interpreted and influenced by X.

The interaction detector is used to identify the interactions between two independent variables (X1 and X2) and the dependent variable Y. q(X1 ∩ X2) is calculated and compared with q(X1) and q(X2) to evaluate different relationships (Table 3). There are five possible interaction relationships: nonlinear-weaken, univariate-weaken, independent, bivariate-enhance, and nonlinear-enhance.

The geographical detector needs categorical variables as the input data. It is necessary to discretize numerical variables into several strata. It is also important to select discretization methods and divide the number of classifications [36]. Having compared the discretization methods, including the equal interval method, the quantile classification method, and the natural breaks method, the latter was selected to classify these variables into 11 categories as the optimal scheme based on the maximum q value.

3.1 Relationships between the RC and driving factors at the watershed scale

The relationships between RC and the 13 driving factors of 29 watersheds are shown in Figure 3. The largest coefficient of determination (R ²) was observed between RC and P (R ² = 0.826, p = 0.000), which showed a positive trend. The lowest coefficient of determination was found between RC and T (R ² = 0.049). Both PET and AET greatly influenced on RC. Slope and TWI rather than elevation were the important topographical factors affecting RC. Of 4 soil factors, only silt was significantly correlated with RC (p < 0.01). The proportions of the 3 main land use types were significantly correlated with RC (p < 0.01). Furthermore, RC was more significantly affected by NDVI (R ² = 0.41, p < 0.01). At the watershed scale, climate was found to predict RC well, followed by topography, land use type, NDVI, and soil.

Figure 3

Scatter plots of runoff coefficient values against precipitation (a), temperature (b), AET (c), PET (d), elevation (e), slope (f), TWI (g), sand (h), clay (i), silt (j), K _S (k), NDVI (l), forest (m), cropland (n), and construction land (o) of Shayin River (SYR) (watersheds are indicated by red dots), Hongru River (HRR) (green dots), Huaihe River mainstream (HHRM) (blue dots), and Shipi River (SPR) (yellow dots). The black dashed lines represent the linear regression of all plots.

3.2 Estimating the mean annual RC map

The mean annual precipitation map and the mean annual AET map were used as inputs in the water balance equation to estimate the mean annual RC map of the study area (Figure 4). RC decreased progressively from south to north. RC in the mountain area was higher than in the plain. The estimated RC values of 29 watersheds were extracted (no data pixels were not considered). Figure 5 presents an exponential relationship between estimated RC and actual RC for 29 watersheds, which yielded R ² and root mean square error (RMSE) of 0.656 and 0.007, respectively. Subsequently, the exponential equation was used to correct the estimated RC map.

Figure 4

The spatial distribution of mean annual runoff coefficient estimated by the water balance equation across the upper reaches of Huaihe River (no data pixels were covered by surfaces of water body, urban and wetland).

Figure 5

Plots of the estimated and the actual runoff coefficients for Shayin River (SYR) (watersheds are indicated by red dots), Hongru River (HRR) (green dots), Huaihe River mainstream (HHRM) (blue dots), and Shipi River (SPR) (yellow dots). The black dashed line and the black solid line represent the linear and exponential regressions of all plots. R ²: coefficient of determination; RMSE: root mean square error.

3.3 Relative influences of driving factors on the RC at the pixel scale

4.1 The RC map

The water balance equation was used to produce a RC map using the mean annual precipitation and the mean annual AET. The mean annual precipitation map came from the spatial interpolation technique using meteorological observation data. The mean annual AET map came from Moderate resolution imaging spectroradiometer (MODIS) global ET products which reached 86% accuracy evaluated by global flux measurement data [11]. Several studies have shown that MODIS ET data could be used for calibrating biophysical models [37], analyzing ET at different temporal and spatial scales [38], computing monthly water balance components [39], and estimating mean runoff over the years [15]. In this study, the mean estimation RC of 29 watersheds was 0.39, close to the actual RC of 0.37. A fairly significant exponential relationship was observed between estimated RC and actual RC from 29 watersheds (Figure 5). The estimated RC values in the Shayin River were higher than the actual RC values, while the estimated RC values in the Shipi River were lower. The main reason is that the MOD16 ET was overestimated from north to south in Chinese river basins [40]. The uncertainties of MODIS ET data were mostly seen in forests and mountains [41]. It is possible that the MODIS ET product did not take into account the topographical relief condition because of the algorithm limitations and uncertainties from input data [11]. It should be noted that RC values in the regions of water body, wetland, and urban and built-up areas could not be estimated. The main reason is that MOD16A3 ET products used MODIS Leaf area index (LAI) and Fraction of photosynthetically active radiation (FPAR) products (MOD15) as the input data. Non-vegetated pixels in these regions could not calculate LAI/FPAR.

4.2 Driving forces of the RC

Climate and watershed characteristics both control the hydrological process of a watershed and determine the RC. The main driving factors of the RCs could be related to scale. The geographical detector is a spatial statistical method to assess the relationships of different geographical strata and reveal the driving forces behind spatial variations [42,43,44].

Climate is the main factor controlling the water cycle [22,45]. The spatial variabilities of climate factors caused significant changes in our study. Precipitation, temperature, and AET decreased from south to north. The RC values of the 29 watersheds also declined from south to north in the study area. High precipitation makes increased runoff more likely. Precipitation proved to be the most important variable for the RC at the watershed scale and across the whole study area. By contrast, at the scale of a single river basin or smaller, the dominant factor which exhibited the highest explanation rate was AET. The main reason is that the spatial variation of precipitation in a small region is relatively slight. In the temperate and subtropical monsoon climate zones, precipitation and evapotranspiration were the main factors affecting runoff [20].

Topography also exerts a strong influence on hydrological processes [25]. Of the topographical variables considered in this study, slope and TWI were more important than elevation at the watershed scale mainly because they could affect the movement of water across watersheds. Mean basin slope could be used to predict RC in various environments [46]. TWI is used to calculate topographic control on wetness and runoff. Runoff generation can be modeled using TWI [47,48]. By contrast, from the watershed scale to the pixel scale, elevation was the most important topographical variable.

NDVI can be used as an important factor to regulate water resources both at the watershed scale and at the pixel scale. NDVI had higher q values. NDVI is a comprehensive reflection of climate, topography, land use, and soil [20]. Higher NDVI values were mainly distributed in forest and mountain regions. The forest density variable was one of the most important factors in predicting the RC [23].

The RC was also found to be highly influenced by land use at the watershed scale. The RC generally increased with an increase in the forest coverage. On the one hand, forest areas can hold more water and produce more runoff than other land use types [30,49]. On the other hand, in this study area, the forest is generally located in mountainous areas, where precipitation is more than 1,024 mm. Cropland coverage was negatively associated with the RC, because of low elevation and precipitation interception. In addition, the irrigation of croplands in the study area primarily uses water from the river systems [49]. Construction land was found to have a negative effect on runoff generation at the watershed scale (Figure 3o), which was not consistent with previous studies [50,51,52]. Some studies found that there is a threshold effect of the impervious area on the RC. When the percentage of construction land area was more than 10% [53] or 20% [50], the construction land area would have an impact on the RC. This study’s largest percentage of construction land areas is less than 18%. The proportion of forest had a significant negative correlation with the construction land. Therefore, forest coverage was the most important land cover type on the RC.

However, land use types were found to have less effect on the RC at the pixel scale. The q values of land use were low. Soil factors and the land use type had only a minimal influence on RC distribution. Soil properties did not seem to exert a major influence on RC. Only silt content was highly correlated with RC (p < 0.01). Moreover, the q values of soil properties were low. Land use, soil types, and geology had a minor control on RCs of the watersheds [24].

4.3 Analysis of interaction between factors

Given that the interactions between different driving factors and the RC are complicated, the interaction detector was used to analyze the interactions between any two driving factors. The interaction results indicated that the explanatory powers of any two-factor interaction on the RC was higher than the q values of a single factor, indicating that the spatial heterogeneity of the RC was largely affected by interactions between factors. Climate and watershed characteristics explained the variation in the RC [30].

4.4 Limitations

There are several limitations to the current study. First, remote sensing products can predict RC at the watershed scale, but uncertainties and accuracies of these products should be considered. MOD16 ET products are inevitably subject to errors resulting from input data and retrieval algorithms. Evaporation on different surfaces can create significant variations. Second, RC values in the construction land could not be estimated. It is feasible to estimate the impervious surface area from remote sensing images and simulate surface runoff by the soil conservation service curve number method [52,54]. Third, the interactions among three or more factors should be further explored in the future. Fourth, numerous factors affecting the runoff and those that predominate in different climatic conditions vary. Their interactive influences can be extremely complex.