Characterizing the groundwater storage–discharge relationship of small catchments in China

Surface flow, baseflow, interflow, etc., are the various elements of streamflow throughout the precipitation process. Generally, the streamflow can be divided into direct runoff and baseflow. Direct runoff, such as interflow or surface flow, responds more rapidly to rainfall events, while baseflow, whose primary source is groundwater, responds slowly to rainfall events (Koskelo et al. 2012; Stewart 2015). The significant role in maintaining water quality and water ecological health is played by baseflow. Therefore, investigating baseflow characteristics can assist in estimating groundwater storage variation, groundwater recharge, and water management (Yin et al. 2011; Nalley et al. 2012; Ahiablame et al. 2013; Gao et al. 2015). It can also detect the catchment hydrological processes or improve the streamflow prediction accuracy (Taormina et al. 2015).

However, the lack or poor quality of runoff data in many catchments prevents accurate estimates of α and β. Therefore, a priori parameterized by the recession curves based on measurable catchment characteristics is a solution to estimate the recession parameters.

The baseflow recession process is primarily controlled by climate and underlying surface characteristics (Price 2011). By deriving the analytical solutions from the Boussinesq equation, we see that the soil hydraulic conductivity and its spatial heterogeneity significantly influence baseflow recession (Harman et al. 2009). Topography and its variation also can affect the shape of the recession curve (Fujimoto et al. 2008; Karlsen et al. 2019). In addition, the shapes of the recession curves are strongly affected by the initial soil moisture state of catchments (Biswal & Marani 2014). Factors influencing the baseflow recession process may also affect the recession parameters. The analysis of these factors will help to parameterize α and β.

Some parametric studies have been carried out on the recession process. For example, Ye et al. (2014) concluded that the average slope (θ), mean saturated hydraulic conductivity at the surface (K_s), drainage density (D_d), the vertical decay parameter of the saturated hydraulic conductivity (f), soil porosity (φ), soil depth (d), and aridity index (AI) are essential for the parameterization of recession parameters α and β for 50 catchments of the USA. Mathias et al. (2016) characterized recession parameters of 120 catchments of the UK using the catchment area (A), baseflow index, and potential evapotranspiration through the multiple linear regression (MLR) analysis. Ali et al. (2014) constructed a physically distributed hydrological model to estimate α and β of 50 catchments in the USA. The results showed that the recession parameters could be predicted through catchment slope and soil hydraulic parameters. However, the results of the parameterization, including the equations and factors influencing the recession parameters in these studies, are different. The reason may be that the diversity of climate and underlying surface conditions between the study areas led to various characteristics of baseflow recession (Patnaik et al. 2018). These results cannot be well applied in China. Therefore, it is necessary to carry out a parameterized study on the baseflow recession based on the observed runoff and catchment characteristics data in China to investigate the baseflow.

The main objective of this study is to derive the parameterization of the baseflow recession process from 382 small catchments in China and transplant the parameterized results. The specific goals are (1) to estimate the recession parameters α and β, (2) to analyze the main factors of recession parameters and using the main characteristics of the catchment to parameterize α and β, and (3) to transplant the parametric equations. The research outcomes help clarify the physical mechanism of baseflow recession and provide data support for the rational allocation of water resources in China.

A total of 412 catchments located in Chinese mainland were collected (Figure 1). In this study, 382 catchments with a Q data record of more than 10 years (1968–1988 and 2006–2014) were selected for parametrization, and 30 catchments with a Q record of 1–10 years were used for the transplantation of parameterization. The physical and climatic characteristics of the selected catchments vary widely, making it possible to analyze the baseflow recession process. The A ranges from 34 to 18,211 km², the θ varies from 0.4° to 28.4°, and the AI varies from 0.11 to 9.43.

The daily Q data were obtained from the ‘Hydrological Yearbook of the People's Republic of China’ which is officially published. The parameters of soil hydraulic characteristics were derived from the China Soil Hydrological Dataset (CSHD) (http://westdc.westgis.ac.cn) (Dai et al. 2013). The soil texture data were taken from the dataset provided by the Data Center for Resources and Environmental Sciences, Chinese Academy of Sciences (RESDC) (http://www.resdc.cn/). The CSHD includes 12 parameters such as field capacity (FC) and saturated hydraulic conductivity. The vertical variation of soil hydraulic characteristics was captured by seven layers to the depth of 1.38 m (i.e., 0–0.05, 0.05–0.09, 0.09–0.17, 0.17–0.29, 0.29–0.49, 0.49–0.83, and 0.83–1.38 m) (Dai et al. 2013). The hydraulic parameters are calculated by the Pedotransfer Functions (PTFs). The inputs of PTFs consist of characteristic parameters of the soil, such as the organic matter content. The hydraulic characteristics data adopt the NetCDF format with a resolution of 1 km × 1 km, which can better capture the spatial heterogeneity of hydraulic parameters (Dai et al. 2013). The soil texture dataset compiled at a scale of 1:1,000,000 consists of three types: sand, clay, and silt data.

The climate, vegetation, topography, and land use data were taken from the RESDC. Based on monthly rainfall and mean temperature data from 1915 stations in China, the climate data such as AI with a spatial resolution of 500 m × 500 m were derived by interpolation method. The normalized difference vegetation index (NDVI) spatial distribution data in 2000 with a resolution of 1 km × 1 km were used as vegetation data. The remote sensing data in 2000 with a resolution of 1 km × 1 km were used for extracting land use information [proportion of water area (PW), proportion of cultivated land area (PC), proportion of woodland area (PF), proportion of grassland area (PG), and proportion of residential land area (PCS)]. The digital elevation model (DEM) was used to extract the terrain data [A, θ, and coefficient of variation of slope (CVS)].

The recession periods should be extracted primarily to estimate the recession parameters, and then the recession curve is used to fit the recession data linearly.

The main criteria for extracting the recessions from the streamflow data are the recessions are extracted from the monotonically decreasing parts of the 2-day moving average of a hydrograph. These recessions must have a minimum length of 4 continuous days. The decline in Q for two consecutive data values has to be smaller than 30% to eliminate the interference from surface runoff (Vogel & Kroll 1992; Ye et al. 2014). The purpose of moving average is to reduce the influence of horizontal scattering on estimating recession parameters caused by low dQ/dt values (Stoelzle et al. 2013). Figure 2 shows an example of the recession periods identified by the above criteria.

The values of catchment characteristics are extracted to parameterize α and β. The specific steps are as follows:

• DEM has been used to extract the A and θ values from the catchment boundaries.

• The catchment boundaries were used to crop the TIFF (Tag Image File Format) images of climatic and physical characteristics, then the catchment characteristic data were calculated. First, for the convenience of calculation, the hydraulic characteristic data in the NetCDF format were transformed into TIFF images. Then, the TIFF images of AI, NDVI, land use, soil hydraulic, and soil texture were cropped by catchment boundaries. Finally, the average and standard deviation of these characteristics in each catchment were calculated.

• Calculating CVS and f values of each catchment. The CVS is the ratio of slope standard deviation to the θ. The calculation of exponential decay parameter, f, refers to Ye et al. (2014), which assumed that the saturated hydraulic conductivity decreases exponentially with depth:

  

  (8)

  where K_d is the saturated hydraulic conductivity at depth d.

The Spearman correlation method is a measure of correlation between variables. This method has been used in many hydrological research works, such as calculating correlation coefficients between the catchment characteristics (Ali et al. 2012). In this study, the Spearman correlation analysis was performed on the above 26 catchment characteristics to select characteristics that are not highly correlated with each other. The details are as follows: (1) calculating the correlation coefficient R between characteristics and recession parameters; (2) for characteristics that are highly correlated with each other (|R| > 0.7), select the ones that have the highest absolute correlation coefficient with recession parameters (Mathias et al. 2016). After the above steps, 12 characteristics (|R| < 0.7) are finally selected to parameterize the recession process. They are A, CVS, FC, AI, PF, PW, PC, PCS, K_s, f, NDVI, and PG.

We selected 287 catchments as training samples, and the remaining 95 catchments were used as test samples to verify the results. The stratified random sampling method was used to divide training and test samples to ensure the uniform distribution and representativeness of samples throughout the population (Cochran 1977).

Before parameterization, a natural logarithm processing was performed on the recession parameters and catchment characteristic variables to reduce the influence of dimensionality. Then, the complete subset regression (CSR) was used in the training database to select the main factors influencing the recession parameters. The MLR was used to parameterize α and β. The CSR method tests all potential combinations of variables, and ‘the best’ combination with the largest R² is selected. Then, the parameterized equations of α and β were derived through the MLR between ‘the best’ variable combination and recession parameters.

The larger the R², the better the fitting effects and the greater the degree of explanation for the change in recession parameters. Generally, the acceptable range of R² is 0.5 and above. On the contrary, the smaller the RMSE, the better the fitting effects of the equations (Santhi et al. 2001; Moriasi et al. 2007; Ahiablame et al. 2013).

The recessions were extracted from the monotonic decreasing parts of a hydrograph based on the criteria mentioned above. The recession parameters α and β were estimated from the recession data by the recession curve fitting linearly. The recession parameters α and β of 382 catchments in China are shown in Figure 4. The recession parameter α range from 3.07e-05 to 0.70, with an average value of 0.06 (Figure 4(a)). Approximately 99% of the α values are between 0 and 0.20. In general, the higher the value of α, the smaller the number of catchments.

The statistical distribution pattern of β values is quite different from that of α values (Figure 4(b)). The values of β range between 0.57 and 3, with an average value of 1.39. Approximately 82% of β values are between 1 and 2. In general, the number of catchments increases first and then decreases with β value increasing. The number of catchments is the largest when the β value is about 1.5.

The spatial distribution patterns of recession parameters are also analyzed to explore the influencing factors of α and β. As shown in Figure 5(a), the α values of the humid climate of southern China are generally lower (α ≤ 0.07) than that of the semi-humid or semi-arid climate of northern China (α > 0.07). The distribution pattern of β values is opposite to that of α values (Figure 5(b)). The β values are greater in humid catchments (β > 1.27) but smaller in semi-humid and semi-arid catchments (β ≤ 1.27). In general, the α values are greater in the arid catchments, while the β values are greater in the humid catchments. In the special case of β = 1, the groundwater outflow behaves as a linear reservoir. By analyzing the spatial distribution of β values and the hydrogeological map, it is found that the main type of aquifer of the catchments with β approximately 1 is the clastic fissure aquifer.

In addition to topography and climate, the other important factors affecting the recession parameters are soil, vegetation, and land use. Therefore, the correlations between the 12 independent characteristics and recession parameters were analyzed. The results show that α and β are significantly correlated with CVS, FC, PC, PF, PW, and PCS (p < 0.05) but not significantly correlated with A and PG (p > 0.05). In addition, NDVI is significantly correlated with α but not with β, while AI, f, and K_s are correlated significantly with β but not with α.

Parameters α and β were characterized using the MLR to construct equations of recession parameters and catchment characteristics. The results of parameterizations of α and β in the training catchments are shown in Table 1. The validation results of parametric equations in the test catchments are shown in Table 2.

All the training catchments can be classified into four types. The first type includes 185 catchments, representing about 64% of all training catchments. The recession parameters of the first type of catchments can be predicted by the equations α₁ and β₁. The R² values of equations are not less than 0.5, indicating that the fitting effects of the equations α₁ and β₁ are acceptable. The recession parameters of the first type of catchments are mainly affected by A, CVS, and FC. In general, A and FC are significantly negatively correlated with α (p < 0.05), while significantly positively correlated with β. On the contrary, the CVS is significantly positively correlated with α, while significantly negatively correlated with β.

The recession parameters of the second type of catchments can be expressed by equations α₂ and β₂. This type consists of 56 catchments accounted for about 20%. The R² values of equations α₂ and β₂ are 0.56 and showing satisfactory fitting effects. The main factors that influence the recession parameters of the second type of catchments are CVS and FC. The CVS is significantly positively correlated with α and significantly negatively correlated with β, while the FC is significantly negatively correlated with α and significantly positively correlated with β.

Thirty-nine catchments (approximately 14%) are classified as type three, and their recession parameters can be fitted by equations α₃ and β₃. The R² values of equations α₃ and β₃ do not exceed 0.4. The low R² may be due to the presence of other variables in addition to the 12 selected variables that have significant impacts on the baseflow recession processes of these catchments. Due to data limitations, the R² cannot be improved more. However, the relationships between variables and recession parameters are significant, indicating that the equations α₃ and β₃ can express recession parameters of the third type of catchments to some extent. The results show that A, CVS, FC, and AI are major factors of parameter α in the third type of catchments, and A, AI, and PW are important factors of parameter β.

The remaining seven training catchments (accounted for about 2%) could not be parameterized, and no significant relationships were found between recession parameters and variables. The reason may be that the selected factors have limited impacts on the recession parameters of these catchments.

The parametric equations obtained from training catchments were used in other 95 catchments to test the general applicability of these equations. As shown in Table 2, the recession parameters of 47% catchments can be parameterized by equations α₁ and β₁, 32% catchments can be parameterized by equations α₂ and β₂, 17% catchments can be parameterized by equations α₃ and β₃, and 4% catchments cannot be parameterized. The results indicate that the parametric equations can estimate recession parameters for most of the catchments, i.e., the parametric equations have universal applicability.

Furthermore, the spatial distribution of parametric equations (Figure 6) was projected on the hydrogeological map to explore the characteristics of various catchments. It is found that the types of aquifers in different catchments are relatively diverse, including porous aquifers, bedrock fissured aquifers, and karst aquifers. In other words, the aquifer characteristics of each type of catchment are not unique.

Based on the estimated values of the recession parameters α and β, the S–Q relationship parameters a and b can be derived according to Equations (4) and (5). The results show that except for the 20 catchments with β values more than 2, the a and b values range between 0–0.99 and 0.70–106.56, respectively.

The influence of the main variables on the recession parameters is analyzed to understand the main mechanism of the baseflow recession process. Overall, the recession parameter α is significantly negatively correlated with A and FC but is significantly positively correlated with CVS and AI. The recession parameter β is significantly positively correlated with A and FC but is significantly negatively correlated with CVS, AI, and PW. The relationships between recession parameters and main variables can be interpreted by residence time and the velocity of outflow (Mathias et al. 2016).

The greater the values of A and FC, the greater the water storage capacity of the catchment, which generally leads to higher residence times and lower values of α (Mathias et al. 2016). Catchments with higher CVS and AI values have lower residence times, resulting in higher α values (Shaw & Riha 2012).

The relationship between A and FC with parameter β is positively correlated. When A value is large, the recession in stormflow and the recession in baseflow are likely to overlap (Mathias et al. 2016). The baseflow recession curves will become sharp; thus, the β values become higher. The impact mechanism of FC on β is that the higher the value of FC, the greater the water storage capacity of the catchment, which leads to an increase in baseflow recession and makes the β values larger. In addition, the correlations between CVS, PW, AI, and β are negative. The reason may be that the higher the CVS, PW, and AI values, the slower the recession rate, thus making the β values smaller.

Consistent with previous studies, the A, CVS, FC, AI, and PW are essential factors in estimating baseflow recession. They pointed out that the recession parameters are the functions of the catchment topography, soil hydraulic, and climate characteristics (Shaw & Riha 2012; Ali et al. 2014; Patnaik et al. 2018). These factors control the recession process by influencing the storage capacity, the contribution ratio of baseflow on storm–runoff, the recession velocity of baseflow, etc.

The parametric equations were transplanted to predict the recession parameters based on the assumption that hydrological responses are similar if catchment characteristics are similar. The effects of the transplantation are verified by the RE of the fitted values of recession parameters and the correlation coefficient between the fitted and the estimated values of recession parameters.

The validation method of transplantation effects is as follows: first, three sets of parameter equations were used to fit the recession parameters. Then, if the RE of recession parameters values fitted by transplanted equations was the smallest, the transplantation effects were represented by the letter x. If the RE was the largest, the letter y was used to indicate the transplantation effects. If the RE was neither the largest nor the smallest, the transplantation effects were represented by the letter z. Finally, the correlation coefficients between the fitted and the estimated values of recession parameters were calculated to evaluate the fitting effects.

The transplantation results of parametric equations in 30 catchments show that equations can parameterize the recession parameters of 26 catchments. Among these catchments, the catchment numbers with transplantation effect of x, y, and z are 19, 4, and 3, respectively. The correlation coefficients between the fitted and the estimated values of α (Figure 7(a)) and β (Figure 7(b)) are greater than 0.60. In addition, the correlation is significant (p < 0.05). Therefore, the transplantation effects of equations are satisfactory. It shows that the parametric equations can be transplanted to predict the recession parameters of catchments.

This paper presents a parameterized study on the baseflow recession process of small catchments in China. The main conclusions are as follows:

• The recession parameters α and β of small catchments in China are 0–0.70 and 0.57–3, respectively. The values of α and β are mostly between 0 and 0.20 and between 1 and 2, respectively. In addition, the S–Q parameters a and b values can be derived from the relationship between the S–Q parameters and the recession parameters, which provides empirical ranges for the calibration of model parameters.

• At the catchment scale, FC, A, CVS, AI, and PW are the main factors that affect the α and β values. In general, the catchment characteristics related to water storage capacity are significantly negatively correlated with α values and significantly positively correlated with β values. The recession parameters are expressed in equations:

  

• The transplantation effects of the equations in the validation catchments are satisfactory. The correlation coefficients between the fitted values and the estimated values of the recession parameters are greater than 0.6. The transplantation of parametric results is significant for the estimation of model parameters and storage of data-scarce catchments.

The project was supported by the National Key R&D Program of China (2018YFC1508103 and 2018YFC1508102), the National Natural Science Foundation of China (51809173 and 51879136), the Water Conservancy Science and Technology Innovation Project in Guangdong Province (202012), and the Open Research Fund Program of the State Key Laboratory of Hydroscience and Engineering (sklhse-2021-A-01).

The authors declare that the research was conducted without any commercial or financial relationships that could be construed as a potential conflict of interest.

H.Y. (Postgraduate) contributed significantly to data analyses and wrote the manuscript. H.H. (Assistant Research Fellow) designed the study and revised the manuscript. Y.L. (Associate Professor) contributed to the idea and revision of the paper. M.T.(Postgraduate) contributed to the collection and consolidation of data. T.Y. (Postgraduate), Z.W. (Associate Professor), L.C. (Associate Researcher), M.Y.A.K. (Professor), and Z.C. (Postgraduate) contributed to the revision of the paper.

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