Extraction Method of Baseflow Recession Segments Based on Second-Order Derivative of Streamflow and Comparison with Four Conventional Methods

Abstract

Baseflow recession analysis is widely used in hydrological research, water resource planning and management, and watershed hydrogeological research. The first step of baseflow recession analysis is to extract the baseflow recession segments from the hydrograph. Different extraction results lead to different analysis results. At present, the four major recession segment extraction methods applied by hydrologists are mostly based on experience, and there is no clear theoretical basis. Therefore, this study derives a second-order derivation (Sec-D) recession segment extraction method based on the power law relationship between storage and discharge. Moreover, by applying the Sec-D method and the four conventional extraction methods to four hydrological stations in the Tao’er River basin in northeastern China, the differences in the recession segment extraction, determination of basin-wide hydrogeological parameters, and groundwater balance estimation are compared. The results demonstrate that, contrary to the four conventional methods, the Sec-D method can effectively eliminate the early recession stage affected by the surface runoff or rainfall and some streamflow data with more than 1% non-sequential error. The hydraulic conductivity of the four basins estimated by the Sec-D method is between 2.3 × 10⁻⁵–4.9 × 10⁻⁵ m/s, and the aquifer thickness is between 131.2 and 202.5 m. However, the four conventional extraction methods may underestimate (by about 2.5 times) the basin-wide hydraulic conductivity and overestimate (by about 3 times) the aquifer thickness. The groundwater balance elements calculated by the Sec-D method and the four conventional methods present similar intra-annual fluctuation characteristics; the correlation coefficients of daily evapotranspiration calculated by the five methods ranged from 0.7 to 0.95, and those of daily effective groundwater recharge ranged from 0.95 to 0.99. The use of the Sec-D method in baseflow recession analyses is significant for future studies and can be combined with conventional methods.

1. Introduction

Baseflow, also known as low flow or drought flow, represents the streamflow formed mainly by discharge of stored water in a basin (aquifer, soil, lake, swamp, etc.) during the dry season [1,2,3,4]. The shape of a baseflow recession curve is an important feature of a basin. Baseflow recession analysis is widely used in hydrological research, water resource planning and management, and watershed hydrogeological research [5,6]. The main application objectives include the estimation of the long-term groundwater storage trends [1] and watershed groundwater balance [7], baseflow separation [8,9], baseflow regionalization [10], estimation of the seepage of groundwater into leaky sewer systems [11], and determination of basin-wide hydrogeological parameters [12,13]. Extensive applications require an accurate and reasonable recession analysis method as their basis.

Baseflow recession analysis can be understood as a statistical analysis and description of the recession curve for a specific basin based on different storage–discharge relationships [14,15,16]. The recession process of baseflow is influenced by many factors, such as aquifer and vegetation characteristics, evapotranspiration, and human activities [6]. The characteristics of the aquifer determine the groundwater storage capacity of a basin, the characteristics of vegetation affect the intensity of transpiration, and evapotranspiration affects the recession process by directly consuming the groundwater storage. Human activities, such as groundwater exploitation, artificial recharge, and agricultural irrigation, clearly have an impact on natural groundwater storage and change the characteristics of baseflow recession. These influencing factors are present to different degrees in different basins, and show significant fluctuations. These fluctuations ultimately increase the difficulty of baseflow recession analysis and result in the diversity of recession analysis methods.

The process of baseflow recession analysis includes three main stages: (1) Extracting the recession segments (Section 2.1 describes four existing conventional methods). Most of these methods used are based on experience. A method with a clear theoretical basis is derived in this study. (2) Choosing a suitable analytical model. The current widely used theoretical models are: linear reservoir model, $\mathrm{S}=\mathsf{\gamma }\mathrm{Q}$ [17]; nonlinear reservoir model, $\mathrm{S}=\mathsf{\gamma }{\mathrm{Q}}^{\mathsf{\beta }}$ [9]; and power law relationship between $-\mathrm{dQ}/\mathrm{dt}$ and Q, $-\frac{\mathrm{dQ}}{\mathrm{dt}}=\mathrm{a}{\mathrm{Q}}^{\mathrm{b}}$ [2]. These three models are essentially based on the power law relationship between storage (S) and discharge (Q) [18,19]. (3) Determining the optimal model parameters, which can be achieved by three widely used methods: linear regression, lower envelopes, and binning [14]. This study focuses on the first stage, which is the premise of recession analysis, and the following two stages are not covered in detail. The power law relationship model and linear regression parameterization method are used in the current research.

The purpose of extracting the baseflow recession segments is to screen out the recession stage which is not affected by precipitation or surface runoff, and to use it for further recession analysis. Precipitation data can be used to exclude the stages directly affected by rainfall. However, daily precipitation data is not available for all catchments, and the end of rainfall does not necessarily mean the end of surface runoff. Therefore, hydrologists typically identify recession stages initially through the downward portion of the hydrograph (dQ/dt < 0) and then remove some early and abnormal recession stages that may be affected by rainfall and surface runoff. At the same time, they also limit the minimum recession length. However, the four conventional recession extraction methods employ different criteria to remove the stages that may be affected by rainfall and surface runoff (Section 2.1), and they do not all provide a clear theoretical basis for establishing these criteria; that is, the formulation of these criteria is quite empirical. On the other hand, the research by Stoelzle, Stahl and Weiler [14] shows that different recession segment extraction methods will lead to different recession coefficients and corresponding storage. The empirical nature of these conventional methods and the fluctuation of recession analysis results is a source of significant confusion for hydrologists, who cannot be sure of the reliability of each method for recession analysis.

Therefore, based on the characteristics of the second-order derivative (Sec-D) of streamflow and the power-law relationship between dQ/dt and Q, this study proposes a new method for extracting recession segments (Section 2.2). Compared with the four conventional methods, this new Sec-D method has a clear theoretical basis and fewer restrictions. Then, this study analyzes the applicability of the Sec-D method based on streamflow data from four hydrological stations in the Tao’er River basin in Northeast China (Section 3) and compares the results to those from the four conventional methods. The differences between the results of recession segment extraction and the results of recession analysis (Section 4.1), the estimation results of basin-scale hydrogeological parameters (Section 4.2), and the estimation results of basin-scale groundwater balance (Section 4.3) are compared and analyzed. Finally, the robustness of the Sec-D method is briefly discussed (Section 4.4).

2. Theory and Methods

2.1. Conventional Recession Segment Extraction Methods

At present, there are four main extraction methods for the baseflow recession segment, developed by hydrologists (

Table 1. Recession segment extraction methods.

Extraction Methods	Criterion	Minimum Recession Length (days)	Excluding Exterior Parts of Recession Segment	Exclusion of Anomalous Recession Decline
Kir	$\mathrm{dQ}/\mathrm{dt}<0$	5 h	First 3 h	--
Vog	Decreasing 3-d moving average	10	First 30%	$\left({\mathrm{Q}}_{\mathrm{i}}-{\mathrm{Q}}_{\mathrm{i}+1}\right)/{\mathrm{Q}}_{\mathrm{i}+1}>30\%$
Bru	$\mathrm{dQ}/\mathrm{dt}<0$	6–7	First 3–4, last 2	$\mathrm{d}{\mathrm{Q}}_{\mathrm{t}+1}/\mathrm{d}\mathrm{t}>\mathrm{d}{\mathrm{Q}}_{\mathrm{t}}/\mathrm{d}\mathrm{t}$
Aks	$\mathrm{dQ}/\mathrm{dt}<0$	5	First 2	CV > 0.10
Sec-D	$\mathrm{d}\mathrm{Q}/\mathrm{d}\mathrm{t}<0,\left|{\partial }^2\mathrm{Q}/\partial {\mathrm{t}}^2\right|<\mathrm{Z}{}^{\#}$	2	--	--

^h the minimum recession length and excluding exterior parts were used in the MATLAB toolbox HYDRORECESSION; ^# Z is the threshold which can be estimated using Equation (2) based on the values of a, b, and Q₅₀.

): the Kir method [20], the Vog method [21], the Bru method [1], and the Aks method [22]. At an early stage, the Kir method uses all pairs of streamflow data during dry periods (dQ/dt < 0). However, in order to avoid the influence of storms and spurious observations, [23] added some criteria: the recession length should be longer than five days, and three days after the flood peak. The Vog method is based on the assumption that a recession begins when a three-day moving average begins to decrease and ends when a three-day moving average begins to increase. These segments must last longer than 10 continuous days, and the decline in streamflow for two consecutive data values has to be smaller than 30%. Furthermore, the first 30% of every recession segment is excluded to avoid the influence of storm and surface runoff. The Bru method eliminates non-recession parts from the hydrograph using the following criteria: “eliminate all data points with positive and zero values of dQ/dt, and also sudden anomalous ones; eliminate three data points after the last positive and zero dQ/dt, and four data points after major events; eliminate two data points before dQ/dt becomes positive or zero; eliminate data points in a drying sequence, which are suddenly followed by a data point with a larger value of −dQ/dt” [1]. The Aks method assumes that the baseflow recession section satisfies the nonlinear relationship of $\mathrm{S}=\mathsf{\gamma }{\mathrm{Q}}^{0.5}$. This method selects the stage in which the CV (coefficient of variation) value between the simulated streamflow (with the model of $\mathrm{S}=\mathsf{\gamma }{\mathrm{Q}}^{0.5}$) and the observed streamflow is less than 0.10 as the “real” baseflow recession segments. The minimum recession length is five days and the first two days of recession are removed. An overview of the main criteria of each method is shown in Table 1. These four methods use many constraints based on different experiences and assumptions, but none have a clear theoretical basis. This can cause confusion for hydrologists when faced with specific applications.

2.2. The Second-Order Derivative of Streamflow and its Application in Extracting Recession Segment

Nejadhashemi et al. [24] reported that the inflection point (where surface runoff ends and the streamflow is composed entirely of baseflow) on the recession limb can be identified where the second-order derivative (Sec-D) of the streamflow is zero. Moreover, before the inflection point, the second derivative clearly fluctuates; after the inflection point, the volatility is small, and the value tends to zero. The authors applied this method to 12 years of streamflow data from a little-tested watershed in the United States. The results indicated that the method is accurate about 87% of the time, which represents a significantly higher success rate than that found in the results of the method developed by Wittenberg [9].

To more intuitively observe the characteristics of the second-order derivatives (Sec-D) of the streamflow, Figure 1 shows a schematic diagram of the Sec-D. The Sec-D of the streamflow presents clear fluctuation characteristics during the flood peak period, but presents little or no fluctuation during the dry season. From a mathematical perspective, the Sec-D of the streamflow reflects the change in the rate of recession or rise. During a flood event, the streamflow is affected by precipitation and surface runoff, and thus, the Sec-D presents obvious fluctuations. However, during recession, the streamflow mainly rises due to the discharge of the aquifer (ignoring the influence of human activities) and satisfies the power law relationship between storage and discharge. Therefore, the Sec-D presents little fluctuation and tends to zero. Based on these characteristics, in this paper we attempt to develop a method for extracting the baseflow recession segment based on the Sec-D of the streamflow.

Brutsaert and Nieber [2] reported the power law relationship between dQ/dt and Q,

$$\frac{\mathrm{dQ}}{\mathrm{dt}}=-\mathrm{a}{\mathrm{Q}}^{\mathrm{b}},$$

(1)

where Q is the baseflow (L³T⁻¹), t is the time (T), a and b are constants.

The derivative of Equation (1) can be obtained as follows:

$$\frac{{\partial }^2\mathrm{Q}}{\partial {\mathrm{t}}^2}=-\mathrm{a}\mathrm{b}{\mathrm{Q}}^{\mathrm{b}-1}\frac{\mathrm{dQ}}{\mathrm{dt}}={\mathrm{a}}^2\mathrm{b}{\mathrm{Q}}^{2\mathrm{b}-1}$$

(2)

For a specific basin in a natural state, a and b are usually regarded as fixed values. Previous studies have demonstrated that b is usually equal to 1 or between 1 and 3 [2,3,5,8,9]. Thus, Equation (2) shows that the Sec-D is positively correlated with a streamflow in the stage of recession. If a, b, and the maximum baseflow (Q_max) are known, the maximum of the Sec-D in the baseflow recession process can be determined, and the part of the Sec-D that is below the maximum value (which can be seen as a threshold Z) can be selected as the baseflow recession segment. The threshold Z here reflects the maximum value that Sec-D can reach during the baseflow recession process. When the Sec-D is greater than this threshold, it indicates that the runoff process is also affected by factors other than baseflow recession. When the Sec-D is less than this threshold, it can be considered that the runoff process is only affected by the baseflow recession process. According to these characteristics, a baseflow recession segment extraction method based on the Sec-D can be obtained (Table 1). It is noted that for more than two consecutive days, the first derivative is less than zero, and the absolute value of the Sec-D is less than the threshold Z.

In order to determine the threshold Z, we need to determine three parameter values: a, b, and Q_max, and then substitute them into Equation (2). The parameters a and b can be determined by conventional recession analysis methods. Previous studies demonstrated that when the streamflow is below the median flow (Q₅₀), it is mainly controlled by the baseflow [6,16]. Thus, herein, a Q₅₀ approximation is used instead of Q_max. The first and second derivatives of the streamflow can be approximated by Equations (3) and (4):

$$\frac{\mathrm{dQ}}{\mathrm{dt}}\approx \frac{{\mathrm{Q}}_{\mathrm{i}}-{\mathrm{Q}}_{\mathrm{i}-1}}{\Delta \mathrm{t}}$$

(3)

$$\frac{{\partial }^2\mathrm{Q}}{\partial {\mathrm{t}}^2}\approx \frac{{\mathrm{Q}}_{\mathrm{i}+1}-2{\mathrm{Q}}_{\mathrm{i}}+{\mathrm{Q}}_{\mathrm{i}-1}}{{\left(\Delta \mathrm{t}\right)}^2}$$

(4)

This method has a clear theoretical basis and few restrictions. The main uncertainty that may cause confusion is the use of Q₅₀ to estimate the Q_max. Another concern is that small measurement errors in the streamflow may significantly affect the calculation results of the second-order derivative. Section 4.4 analyzes the rationality of Q₅₀ and the impact of measurement errors in detail.

2.3. Determination of Basin-Wide Hydrogeological Parameters Based on Recession Analysis

Brutsaert and Nieber [2] derived three analytical solutions of the one-dimensional Boussinesq equation which describes the drainage of water from an ideal aquifer (a non-confined, rectangular prism with a characteristic breadth B) under Dupuit’s assumption. The solutions have the same form as Equation (1): $\frac{\mathrm{dQ}}{\mathrm{dt}}=-\frac{1.1337}{{\mathrm{K}\mathsf{\mu }\mathrm{D}}^3{\mathrm{L}}^2}{\mathrm{Q}}^3$ for the short term, $\frac{\mathrm{dQ}}{\mathrm{dt}}=-\frac{4.8038{\mathrm{K}}^{0.5}\mathrm{L}}{{\mathsf{\mu }\mathrm{A}}^{1.5}}{\mathrm{Q}}^{1.5}$ for the long term, and $\frac{\mathrm{dQ}}{\mathrm{dt}}=-\frac{0.3465{\mathsf{\pi }}^2{\mathrm{KDL}}^2}{{\mathsf{\mu }\mathrm{A}}^2}\mathrm{Q}$ for both short and long terms. The constant a in Equation (1) is directly proportional to the basin-wide groundwater parameters: for the short term, a₁ = $\frac{1.1337}{{\mathrm{K}\mathsf{\mu }\mathrm{D}}^3{\mathrm{L}}^2}$, b₁ = 3; for the long term, a₂ = $\frac{4.8038{\mathrm{K}}^{0.5}\mathrm{L}}{{\mathsf{\mu }\mathrm{A}}^{1.5}}$, b₂ = 1.5. Hydrogeologists can further calculate the basin-wide hydrogeological parameters by determining the values of a₁ and a₂ based on the observed relationship between dQ/dt and Q [2,25,26]. Subsequently, Parlange et al. [27] presented a single analytical formulation that yields a smooth transition between the short term flow regime, with b = 3; and the long term flows, with b = 1.5. This formulation is based on the non-dimensional definitions for discharge, ${\mathrm{Q}}^{*}=\frac{{\mathrm{B}}^2}{{\mathrm{KAD}}^2}\mathrm{Q}$, and time, ${\mathrm{t}}^{*}=\frac{\mathrm{KD}}{{\mathsf{\mu }\mathrm{B}}^2}\mathrm{t}$. Parlange et al. [27] also found that the transition point between the short- and long term flow regimes is at t* = 0.5625, and the coordinates of the dimensionless transition point are log Q* = −0.1965, log (|dQ*/dt|) = 0.0918. Then, they related the translation magnitudes in the horizontal or abscissa (log(Q)) and vertical or ordinate (log(|dQ/dt|)) directions between the observed transition point and the dimensionless transition point, i.e., H and V, respectively, to the basin-wide hydrogeological parameters as follows:

$$\mathrm{K}={\mathrm{V}\mathsf{\mu }\mathrm{B}}^2/\left(\mathrm{HD}\right)$$

(5)

$$\mathrm{D}={\mathrm{H}}^2/\left(\mathrm{VA}\mathsf{\mu }\right)$$

(6)

$$\mathrm{B}=\mathrm{A}/2\mathrm{L}$$

(7)

where K is the hydraulic conductivity, D is the aquifer thickness, μ is the drainable porosity, B is the aquifer breadth, A is the basin area, and L is the total channel length. Mendoza, Steenhuis, Walter and Parlange [12] outlined the development process of the above theories and methods in detail. The drainable porosity, μ, generally exhibits a smaller, more predictable range of values than the other parameters, such as between 0.00001 and 0.001 for fractured rock [12]. In practical applications, the value or range of μ can be given according to the characteristics (mainly the lithology and fissure development degree) of the actual watershed. B, A, and L can also be determined at the same time, and then the values or ranges of K and D can be calculated according to Equations (5) and (6). Oyarzún, Godoy, Núñez, Fairley, Oyarzún, Maturana and Freixas [13] reported that the intersection of the lower 10% regression line with a slope of 3 and the least squares regression line for log (|dQ/dt|) as a function of log Q data can reflect the average position of the transition point. Herein, this approach will be followed.

2.4. Groundwater Balance Estimation Based on Recession Analysis

Based on the nonlinear reservoir model, Wittenberg and Sivapalan [7] presented a process for estimating the basin-scale groundwater balance based on baseflow separation and recession analysis. In this process, the groundwater discharge (Q) is first determined through the separation of baseflow, then the evapotranspiration (ET) is estimated based on the deviation of monthly recession coefficients and the baseflow sequence, and finally the effective groundwater recharge (GWR) is estimated by the groundwater balance equation based on the continuous baseflow and evapotranspiration sequences. For a detailed introduction, please refer to Wittenberg and Sivapalan [7]. This study does not repeat the discussion; however, the equations used in the three main stages of the process are described in Appendix A.

3. Study Area and Application

3.1. Study Area

The Tao’er River basin is located in northeastern China, with an elevation ranging from 100 to 1600 m above sea level (Figure 2). The Tao’er River originates from Gaoyue Mountain (which is part of the Greater Khingan Mountains) and flows from northwest to southeast with a total channel length of 563 km. The topographic conditions are variable and the slopes of the valleys decrease from upstream to downstream, varying from 1.66% to 0.02%. The drainage area cover 41,200 km², of which the bed rock mountain area accounts for approximately 65%, and is mainly distributed above the Zhenxi Hydrological Station [28,29]. Coverage of forest and grassland in the basin is approximately 25% and 36%, respectively, and the remainder of the area is mostly cultivated with corn [29]. The Tao’er River basin is in a temperate continental monsoon climate zone, with an average annual rainfall of 463.6 mm, and an average annual temperature of −2.7 to 5 °C. Precipitation varies by month, with July and August accounting for approximately 50% of the annual rainfall [30].

According to existing geological and hydrogeological survey results, it is known that the Tao’er River basin study area is widely distributed with Mesozoic and Cenozoic volcanic rocks, volcanic clastic rocks, and granite. Bedrock fissure aquifers are distributed across most of the basin, and pore phreatic aquifers are distributed only in river valleys. The pumping test results of Jia [31] in the valley area to the south of Wulanhaote city show that the drainable porosity (μ) is between 0.00027 and 0.0022. This range is the result of the mixing action of the valley pore phreatic aquifer and the bedrock fissure aquifer. The drainable porosity (μ) of the pore phreatic aquifer is usually large, so the μ of the bedrock fissure aquifers widely distributed outside the valley will be lower than these values. This also shows that the range of μ in this research basin is largely consistent with the range (0.00001–0.001) recommended by Mendoza, Steenhuis, Walter and Parlange [12].

This study is based on the daily streamflow observation data of Suolon, Chaersen, Dashizhai, and Zhenxi hydrological stations (Figure 2) in the basin from 1964 to 1989. Subsequent baseflow recession segment extraction, recession analysis, and estimation of hydrogeological parameters and basin-scale groundwater balance were all based on these streamflow data. The unpublished streamflow data were provided by Songliao Water Resources Commission, Ministry of Water Resources. The characteristics of the drainage area and length of the channel controlled by the four hydrological stations are shown in

Table 2. Basic characteristics of the four hydrological stations.

Hydrological Station	Available Daily Streamflow Data	A	L	B	Q50	Z #
		km2	km	km	m3/s	m3/(s⋅d2)
Chaersen	1964–1989	7872	285	13.8	8.44	0.015
Dashizhai	1964–1989	7656	215	17.7	5.90	0.011
Suolun	1964–1989	5893	230	12.8	7.91	0.017
Zhenxi	1964–1989	18,462	628	14.7	10.70	0.020

^# Z is estimated using Equation (2) based on the values of a (Kir), b (Kir), and Q₅₀.

.

3.2. Application

Based on the daily streamflow observation data of the four hydrological stations from 1964 to 1989, the baseflow recession segments were screened by the five methods shown in Table 1. The differences between the extraction results of different methods in the same hydrological station are compared and analyzed in Section 4.1. The estimation of the threshold (Z) of the Sec-D method is shown in Table 2.

The log (|dQ/dt|) vs. log Q diagrams are drawn according to the extraction results from 1964 to 1989, and the lower 10% regression line with a slope of 3 and the least squares regression line are presented. Hydrogeological parameters were calculated using the method described in Section 2.3, and the results of the hydraulic conductivity (K) and aquifer thickness (D) under different extraction methods are compared and analyzed in Section 4.2. As mentioned earlier, the range of drainable porosity (μ) in this study is basically consistent with the recommended range of Mendoza, Steenhuis, Walter and Parlange [12]. Therefore, this range (0.00001–0.001) was still used in this study but, for comparison purposes, this study only calculated the hydrogeological parameters when the drainable porosity (μ) is equal to the range average 0.0001.

Based on the recession analysis results ($\mathsf{\gamma }=1/\left(2\mathrm{a}-\mathrm{a}\mathrm{b}\right)$ and ($\mathsf{\beta }=2-\mathrm{b}$) of multi-year (1964–1989) monthly average and multi-year (1964–1989) annual average, the groundwater balance of the four basins controlled by the four hydrological stations was calculated by the method described in Section 2.4. The groundwater balance results are expressed in the form of a multi-year daily average, as reported by Miller et al. [32], and the differences between different extraction methods are compared and analyzed in Section 4.3. In order to facilitate the comparative analysis, the β value of each hydrological station was fixed as the average of the five extraction methods (Chaersen 0.76, Dashizhai 1.09, Suolun 0.86, Zhenxi 0.91). The monthly recession analysis of the four conventional extraction methods was completed with the MATLAB toolbox HYDRORECESSION (Arciniega-Esparza et al., 2017). The Sec-D of the streamflow was approximated by Equation (4), and the monthly recession analysis of the Sec-D extraction method was completed using Origin (https://www.originlab.com/, Northampton, MA, USA).

4. Results and Discussion

4.1. Comparison of the Recession Extraction and Analysis

The recession segments extracted by the five recession extraction methods were statistically calculated, and the Kir method was used as the reference to calculate the extraction proportion. Because the Kir method has the most relaxed restrictions, it essentially retains all of the recession stages. The results are shown in

Table 3. Extraction proportion of different recession extraction methods (Kir method as a reference).

Hydrological Station	Extraction Proportion
	Kir	Vog	Bru	Aks	Sec-D
Chaersen	100%	87%	68%	66%	61%
Dashizhai	100%	79%	63%	60%	62%
Suolun	100%	92%	69%	63%	64%
Zhenxi	100%	85%	68%	60%	56%

. For the same hydrological station, the Vog method presents the largest extraction proportion, followed by the Bru method, and the Aks or Sec-D methods have the smallest proportion. Using the Kir method as a reference, the diagrams of −dQ/dt as a function of Q of the recession segments extracted by the four methods are plotted (Figure 3, Figures S1–S3). The results demonstrate that there is no significant difference in the distribution characteristics of the extracted segments of the Kir, Vog, Bru, and Aks methods, while the segments extracted by the Sec-D method are concentrated in the stage of low streamflow. To observe the extraction results more clearly, the segments extracted by the different methods were demonstrated by taking the recession stage of Suolun Station from August to December 1985 as an example (Figure 4). Figure 4 shows that the segments extracted by the Sec-D method clearly exclude many early stages of recession, but the segments extracted by the Aks method clearly exclude some late stages of recession. The Aks method’s hypothesis of $\mathrm{S}={\mathsf{\gamma }\mathrm{Q}}^{0.5}$ may not be applicable to the basins of the current study because the exponents of these basins are close to 1 (

Table 4. Recession analysis results of the four hydrological stations under different extraction methods.

Hydrological Station	a					b
	Kir	Vog	Bru	Aks	Sec-D	Kir	Vog	Bru	Aks	Sec-D
Chaersen	0.024	0.022	0.022	0.024	0.022	1.23	1.25	1.25	1.22	1.23
Dashizhai	0.053	0.052	0.053	0.048	0.051	0.91	0.94	0.87	0.97	0.85
Suolun	0.033	0.031	0.033	0.031	0.032	1.14	1.16	1.13	1.15	1.13
Zhenxi	0.033	0.027	0.029	0.030	0.034	1.09	1.12	1.11	1.11	1.01

, β = 2 − b). No significant difference between Bru, Vog, and Kir method was found. Thus, the Sec-D method eliminates the stage that is clearly affected by surface runoff or rainfall in the early stage of the recession, whereas the four conventional methods do not effectively eliminate them. Although the extraction ratios of the Bru and Aks methods are smaller, the removal effect of the stages affected by surface runoff is not obvious.

The recession coefficients (a, b) under different extraction methods were calculated by the least squares method. The a values of the four hydrological stations ranged from 0.022 to 0.053, and the b values ranged from 0.85 to 1.25 with an average of 1.09 (Table 4). The recession coefficients of the different extraction methods at the same hydrological station do not exhibit significant differences. This shows that although the Sec-D method eliminates more early recession stages, it does not have a significant impact on the multi-year recession analysis results (recession coefficients a and b). Furthermore, it also does not significantly affect the application results based on these coefficients, such as estimating basin-scale groundwater storage [1], reverse baseflow separation (Appendix A1), and baseflow regionalization [10]. However, Stoelzle, Stahl and Weiler [14] observed that there are significant differences in the recession characteristics obtained by different extraction and recession analysis methods. The reasons for the small differences may be the following: First, only the least squares method is used to calculate the parameters, without comparison to other methods. Second, the basins studied in this paper are in a semi-arid area, while the basins studied in Stoelzle, Stahl and Weiler [14] are in a humid area. In semi-arid basins, the rainfall is relatively low, and the baseflow recession segments affected by rainfall and surface runoff are relatively short, which may lead to similar extraction effects to those of different extraction methods. Third, this study is based on four nested basins, without comparison to other regions and climate conditions.

4.2. Comparison of the Results of Estimating Basin-Wide Hydrogeological Parameters

According to the method described in Section 2.3, the hydrogeological parameters (K, hydraulic conductivity and D, aquifer thickness) under different extraction methods were calculated, and the drainable porosity was fixed at 0.0001. The results are shown in Table S1 and Figure 5. As mentioned in Section 4.1, there is no significant difference between the extraction results of the four conventional extraction methods at the same hydrological station. Therefore, there is no significant difference at the transition point of the baseflow recession. With the exception of the obvious fluctuation of K at Dashizhai station, there is no significant difference in the hydrogeological parameters calculated by the four conventional methods at the same hydrological station. However, K calculated by the Sec-D method is larger (by about 2.5 times) than that calculated with the four conventional methods, and D is smaller (by about 3 times) than that calculated by the four conventional methods (Figure 5).

The key to the estimation method of the basin-wide K and D described in Section 2.3 is the determination of the transition point of log (|dQ/dt|) as a function of log Q. The position of the transition point will directly affect the estimation results of K and D. The positions of the transition points of the four conventional methods in the same basin are effectively the same (Figure 3, Figures S1–S3), so the estimation results of K and D are similar. However, the transition point of the Sec-D method is clearly shifted to the low flow side, and the estimation results of K and D are also significantly different from the four conventional methods. This shift is due to the fact that the Sec-D method removes the early stage of the recession, which may be affected by surface runoff or rainfall.

The pumping test results of Jia [31] in the valley area to the south of Wulanhaote city show that the hydraulic conductivity (K) is between 5.8 × 10⁻⁴ m/s and 5.8 × 10⁻³ m/s. Spitz and Moreno [33] concluded that the empirical value of hydraulic conductivity (K) of weathered granite is between 3.3 × 10⁻⁶ and 5.2 × 10⁻⁵ m/s. The K calculated in this study is the average value at the scale of the overall basin scale and has a theoretical value of between 3.3 × 10⁻⁶ and 5.8 × 10⁻³. From Figure 5, it can be seen that the values of K calculated by the five recession segments extraction methods are all within the above range, and all of these results appear to be reasonable.

Oyarzún, Godoy, Núñez, Fairley, Oyarzún, Maturana and Freixas [13] reported that the mean of the maximum elevation difference (D_a) between the highlands and riverbed on both sides of a river channel reflects the maximum possible value of the aquifer thickness in the basin, and that the real aquifer thickness is close to 0.3D_a in an arid basin in north-central Chile. According to the method developed by Oyarzún, Godoy, Núñez, Fairley, Oyarzún, Maturana and Freixas [13], D_a is estimated by regularly arranging four to six transverse relief profiles in the basin controlled by each hydrological station. The average D_a of each basin is 466 m at Suolun station, 305 m at Dashizhai station, 394 m at Chaersen station, and 362 m at Zhenxi station. It can be observed from Figure 5b and Table S1 that the results of the aquifer thickness based on the four conventional extraction methods are all greater than D_a, that is, the four conventional methods lead to a significantly overestimated aquifer thickness. It can be seen from Equation (5) that the estimated K value is inversely proportional to the D value. If the D value is overestimated, the K value will be underestimated. However, the aquifer thickness based on the Sec-D method is below D_a, and the values of D/D_a are 0.43 in Suolun station, 0.43 in Dashizhai station, 0.44 in Chaersen station, and 0.30 in Zhenxi station. In fact, the thickness D of the aquifer calculated by the method outlined in Section 2.1 reflects the thickness of the rock layer that can be saturated in the early recession [2]. It is usually impossible for the water table to reach this theoretical maximum height D_a, especially in bedrock mountain areas. In addition, with the increase of the depth, the development degree of rock fractures gradually decreases, and it is difficult to ensure that there are connected fractures throughout the thickness of Da in an actual watershed. These ratios obtained by the Sec-D method are slightly larger than that of the arid basin (0.3), which seems reasonable for semi-arid basins. This also suggests that the Sec-D method calculates more reasonable basin-wide hydrogeological parameters than the four conventional methods.

4.3. Comparison of the Results of Groundwater Balance Estimation

The baseflow separation method described in Appendix A1 depends on the recession coefficients $\mathsf{\gamma }=1/\left(2\mathrm{a}-\mathrm{a}\mathrm{b}\right)$ and $\mathsf{\beta }=2-\mathrm{b}$. As mentioned in Section 4.1, the difference between the recession coefficients calculated by the different extraction methods at the same hydrological station is small, and the β value of the same hydrological station is fixed when estimating the groundwater balance. Thus, the difference of the baseflow separation results based on different extraction methods at the same hydrological station is not obvious. The baseflow indices (calculated as the ratio of total baseflow to streamflow, and used in this study to compare the differences in the calculation results of groundwater discharge under different extraction methods) calculated by each of the five extraction methods at the Chaersen, Dashizhai, Suolun, and Zhenxi stations ranged from 0.54 to 0.55, 0.60 to 0.61, 0.52 to 0.53, and 0.53 to 0.56, respectively (Table S2).

The calculation results of the multi-year monthly average recession coefficient (γ) obtained from different extraction methods are shown in Figure 6. The fluctuation trend of γ is similar for each method at the same hydrological station, but there are larger deviations between the Sec-D method and the four conventional methods by comparing the results for different months. In general, γ for each of the four hydrological stations has a similar intra-annual trend. In summer (May to August), the evapotranspiration and agricultural exploitation increase, the baseflow recession is faster, and γ is smaller. In winter, the evapotranspiration decreases, agricultural exploitation stops, baseflow recession slows down, and γ is larger. However, γ for each hydrological station, based on the Sec-D method, presents an unusually large value in one month during the summer, which may be due to the continuous concentrated rainfall in the summer caused by the continental monsoon climate, which significantly reduces evapotranspiration and agricultural exploitation. In addition, γ calculated by the five extraction methods of the four hydrological stations presents an obviously low value in November, which may be due to the rapid decrease of the temperature causing riverbeds to freeze and blocking the drainage of groundwater to the river (the multi-year average temperature in the Tao’er River basin in November is −5 °C).

The method described in Appendix A2 was used to estimate the multi-year daily average evapotranspiration. The amount of evapotranspiration calculated by this method depends on the difference between γ and γ_R. The results are shown in Figure S4. Similar to γ, the fluctuation trend of the evapotranspiration is similar for each method at the same hydrological station, but there are large deviations between the Sec-D method and the four conventional methods when the results for different months are compared. In summer, when γ is significantly larger, the evapotranspiration calculated by the Sec-D method is significantly smaller. In winter, when γ is significantly smaller, the evapotranspiration calculated by the Sec-D method is higher. The average annual evapotranspiration results are shown in Table S2. The results based on the four conventional methods are slightly smaller than those calculated by the Sec-D method.

According to the method outlined in Appendix A3, the multi-year daily average effective groundwater recharge was calculated and is shown in Figure 7. With the exception of the recharge peaks in Suolun and Zhenxi stations in July, the difference in the daily average recharge, based on the five extraction methods, is small. The four hydrological stations present the same recharge mode, which has three distinct recharge stages and an obvious peak value in July. The first recharge stage occurs from March to April, mainly due to the temperature rise and a large amount of ice and snow meltwater infiltration. The second and third stages occur between June and September, mainly caused by concentrated rainfall infiltration due to the monsoon climate. The results of the average annual recharge are shown in Table S2. Compared with the Sec-D method, the annual recharge based on the conventional method is slightly smaller, but the deviation is not obvious.

In order to more clearly observe the difference of daily groundwater balance estimation results of the five methods, we drew the normalized Taylor diagrams of daily evapotranspiration (Figure S5) and effective groundwater recharge (Figure 8) estimated by the five methods (the Sec-D method as the reference series). The normalized Taylor diagram can provide a concise statistical summary of how well patterns match each other in terms of their correlation, their root-mean-square difference, and the ratio of their standard deviation. For a detailed introduction to the Taylor diagram, please refer to [34,35,36]. It can be seen from Figure S5 that the correlation coefficients between the daily evapotranspiration calculated by the four conventional methods and that calculated by the Sec-D method are between 0.7 and 0.95, the root-mean-square differences are between 0.4 and 0.8, and the ratios of standard deviation are between 0.75 and 1.25. Generally speaking, there is a strong correlation between the evapotranspiration estimated by the four conventional methods and that by the Sec-D method, but there are still some deviations. It can be seen from Figure 8 that the correlation coefficients between the daily effective groundwater recharge calculated by the four conventional methods and that calculated by the Sec-D method are between 0.95 and 0.99, the root-mean-square differences are between 0.1 and 0.4, and the ratios of standard deviation are between 0.75 and 1.25. In general, the daily effective groundwater recharge estimated by the four conventional methods has a high correlation with that estimated by the Sec-D method, and the deviations are very small. Thus, the impact on the estimation of basin-scale groundwater balance after removing more early recession stages is not obvious. Furthermore, the results are close to those of the four conventional methods with the exception of large fluctuations in individual months.

4.4. The Robustness of the Sec-D Method

Replacing the maximum baseflow with Q₅₀ (Section 2.2) is largely an empirical consideration. In order to analyze the rationality of Q₅₀, we calculated the threshold (Z) for the different percentages of streamflow (Q₁₀, Q₂₀, …, Q₉₀) of the four hydrological stations. The results are shown in Figure S6. The threshold (Z) of four hydrological stations increases exponentially with the increase of percentage streamflow. When the streamflow is greater than Q₅₀, the Z value of Dashizhai station is significantly different from that of the other three stations, which indicates that there are obvious differences in the rainfall-runoff and recession processes between Dashizhai station and the other three stations. Then, we calculated the extraction proportion under the thresholds of different percentages of streamflow. The results are shown in Figure 9. The extraction proportion of the four hydrological stations shows the same fluctuation trend with the increase of the percentage of streamflow: there is no obvious change in the extraction proportion from Q₁₀ to Q₄₀, which is about 30%; from Q₄₀ to Q₅₀, it increases significantly to about 60%; from Q₅₀ to Q₇₀, it is approximately stable at about 60%; and from Q₇₀ to Q₉₀, it gradually increases to about 80%. Q₁₀, Q₂₀, Q₃₀, and Q₄₀ will clearly underestimate the maximum value of baseflow, making the extraction proportion smaller; Q₈₀ and Q₉₀ will clearly overestimate the maximum baseflow, leading to the extraction results containing more recession sections affected by rainfall runoff. Therefore, Q₅₀ will not clearly underestimate or overestimate the maximum baseflow, and thus represents a reasonable choice. It should be noted that, when Q₆₀ and Q₇₀ are used to replace the maximum baseflow, the extraction results are not significantly different from those using Q₅₀.

The Sec-D method developed in this paper is highly dependent on the calculation results of the second-order derivative of the streamflow. Small errors in the streamflow may significantly affect the second-order derivative and, thus, the segments extracted by the Sec-D method. Two main types of error in streamflow exist: random error and systematic error. To investigate, we analyzed the impact of these two kinds of errors. First, the effect of random errors in a single time step was analyzed. To achieve this purpose, we gradually added random errors between −5% and 5% to the recession segments extracted by the Sec-D method from 1964 to 1989 (using Excel), and then analyzed the relationship between the second-order derivative of the streamflow and the random errors (Figure 10). It can be clearly seen from Figure 10 that, as the absolute value of the random error increases, the |Sec-D| value of the streamflow increases significantly. When the absolute value of the random error is greater than 1%, the |Sec-D| values of the four hydrological stations are significantly greater than their threshold Z (Table 4); that is, when there is a random error greater than 1% in a single time step, the Sec-D method eliminates it. Secondly, the influence of systematic errors in the entire time step was analyzed. To achieve this goal, we added fixed errors of 1%, 5%, and 10% to the recession segments extracted by the Sec-D method from 1964 to 1989, and further analyzed the changes of the second-order derivative. The results show that when a constant proportion of errors is added, the changes of the second-order derivatives are very small, and most remain below the threshold Z (the results are not shown in this paper); that is, the Sec-D method does not eliminate the recession stages with stable systematic errors. It should be noted that these systematic errors do not significantly change the estimation results of dQ/dt. This is because the fixed errors in the two streamflow values cancel each other using Equation (3) to estimate dQ/dt, so do not significantly influence the results of recession analysis.

5. Conclusions

In this study, a Sec-D recession segment extraction method was derived based on the power law relationship between storage and discharge. In addition, the applicability and robustness of the Sec-D method in actual watersheds and the differences from the four conventional methods were analyzed through application of the method to four nested watersheds in northeastern China. The following conclusions can be drawn:

(1)

The Sec-D method has a clear theoretical basis with few restrictions, and maintains an acceptable robustness. Its extraction results can be applied to recession analyses, basin-wide hydrogeological parameter determination, and groundwater balance analyses.

(2)

Among all of the methods investigated, the Sec-D method can effectively eliminate the early recession stage affected by surface runoff or rainfall and streamflow values with more than 1% on-sequential error.

(3)

The elimination effect of the Sec-D method will not have a significant impact on the results of recession analysis (recession coefficients a and b), their direct application (estimation of basin-scale groundwater storage, reverse baseflow separation, and baseflow regionalization), and the estimation results of basin-scale groundwater balance. However, the method will clearly affect the estimation results of basin-scale hydrogeological parameters.

(4)

The Sec-D recession extraction method can be used to calculate credible hydrogeological parameters. The four conventional extraction methods examined may underestimate the basin-wide hydraulic conductivity (K) and overestimate the aquifer thickness (D).

The use of the Sec-D method for baseflow recession analyses is significant for future studies and can be combined with conventional methods, particularly in the estimation of basin-scale hydrogeological parameters. In addition, future research should explore the applicability of the Sec-D method to other countries and climatic conditions.