Evaluation of seasonal catchment dynamic storage components using an analytical streamflow duration curve model

Abstract

Dynamic storage refers to groundwater storage that is sensitive to rainfall infiltration, streamflow generation, evapotranspiration, and other variables involving groundwater gain or loss. It plays a crucial role in habitat maintenance and the mitigation of environmental impacts on regional hydrological behaviors. Dynamic storage can be separated into direct storage, which contributes to the river channel, and indirect storage, which is insensitive to streamflow. The combination of diverse approaches would provide an estimation of the two storage types. This study estimated optimal baseflow coefficients and direct storage in the wet and dry seasons using an analytical streamflow duration curve model in eight catchments of the Choushui River Basin from 2013 to 2017. The water balance approach was then combined to assess indirect storage for evaluating seasonal dynamic storage components. The model applicability for each catchment of the Choushui River Basin in the wet and dry seasons was assessed using the similarity between observed and simulated flow duration curves, namely Kolmogorov–Smirnov distance. We also applied it to assess the performance difference between model and streamflow recession analysis, which is typically used to estimate baseflow coefficients. The results demonstrated that seasonal differences in baseflow coefficients were related to catchment characteristics as well as the aquifer extent through which groundwater flows. The model utilizing maximum likelihood estimation exhibited superior performance than streamflow recession analysis and was highly applicable in our study area in wet and dry seasons. Dynamic storage components demonstrated a considerable difference in the additional groundwater storage between dry and wet seasons and a loss of direct storage was observed in most catchments during the dry season.

1 Introduction

Groundwater resource management and quantification are restricted by monitoring challenges and invisibility of aquifers, causing groundwater overuse in several regions of the world [1]. Groundwater storage is not only a major factor in controlling baseflow physics and chemistry but also plays a role in regulating thickness of the vadose zone, ecological availability of water resources, and evapotranspiration. Previous studies have shown that larger storage capacity results in hydrological connectivity and deep groundwater flow [2, 3]. Rihani et al. [4] mentioned that strong correlations of groundwater table with geomorphology, aquifer heterogeneity, vegetation type, and regional climate demonstrated that seasonal groundwater table changes are controlled by horizontal drainage of the channel and vertical loss from groundwater. Thus, understanding the relationships between groundwater storage, hydrogeological structures, and climate patterns may help improve long-term hydrological predictions under future climate scenarios. Exploring groundwater storage response mechanisms is also vital for comprehending quality and quantity of groundwater, hydrogeological structures, and interactions with ecological or human activities.

Hydrologists have pursued to identify physical attributes that sufficiently describe different basin hydrological behaviors to avoid high parameterization and model uniqueness issues [5]. This has resulted in the development of storage–discharge behavior, which describes the groundwater discharge process from aquifers of catchment or mountainous areas to the river channel [6]. Streamflow recession analysis (REC), proposed by Brutsaert and Nieber [7], has been widely adopted to explore the storage–discharge relationship and to estimate groundwater storage [8,9,10]. Based on the assumption that groundwater discharge dominates streamflow during the long-term dry period, the storage–discharge relationship can be determined by data selection and fitting. Owing to regional climatic conditions and statistical rationality, several different combinations of low-flow selection criteria with different levels of rigor [11, 12] and parameter fitting [13, 14] have been developed utilizing this technique. Jachens et al. [15] suggested that, if recession analysis is integrated with interdisciplinary science, a representative method must still be chosen to establish future development. The analytical flow duration curve (FDC) model developed by Botter et al. [16,17,18] is presumed to be a recession analysis alternative [19]. This model extends the relationship between stochastic rainfall and soil water content to the streamflow generation process, deriving an analytical expression with baseflow coefficients for the daily streamflow probability distribution. The advantage of this method compared with the REC is its availability for all streamflow data for at least one year. It has also been applied in many regions with different hydrological and climatic conditions [20,21,22].

Non-uniqueness of the storage–discharge function makes it challenging to elucidate physical processes, and there may be a huge number of seasonal groundwater storage changes that are unlinked to discharge from aquifers to channel. Conceptual hydrological models also usually characterize these storages as different storage components, which are assumed to be important factors in catchment discharge response and immediate or delayed phenomena of the storage–discharge relationship [23]. Estimating storage components and their temporal changes can offer more insights into catchment hydrological processes and enhance hydrological simulation and prediction. Additionally, due to different definitions of groundwater storage estimated using different approaches, a combination of these methods can help us comprehend the storage components [24]. Dralle et al. [25] proposed that total dynamic storage (S_T) could be separated into two components: direct storage (S_d), which contributes to streamflow, and indirect storage (S_i), which is less sensitive to streamflow. Their results indicate that aquifer properties control the S_T composition. This will help predict the impact of climate change on groundwater evaporation and promote the consideration of groundwater processes in hydrological models.

The aim of this study was to use the analytical FDC model in combination with the water balance method to compute seasonal baseflow coefficients for quantifying dynamic storage components in eight catchments of the Choushui River Basin. Our study also explored model applicability of each catchment in wet and dry seasons using the Kolmogorov–Smirnov distance (c^KS). It represents the maximum distance between the simulated and observed FDC at a specific streamflow and could be appropriate for assessing the similarity between cumulative distribution functions. Although the model is available for all streamflow data, the observed FDC still need enough data to constitute the streamflow distribution characteristics for parameter fitting and model applicability assessment. To address the representativeness of FDC, this study used five years of streamflow data from 2013 to 2017 in each catchment for analysis. We also ascertained the differences in the ratio between direct and indirect storage owing to seasonal rainfall. The results of this study will help enhance our understanding of catchment storage response mechanism and the relationship between seasonal rainfall differences and aquifer properties.

2 Materials and methods

2.1 Study area

The Choushui River Basin, located in the central region of Taiwan (Fig. 1). The Choushui River originates from the main and east peak of the Hehuan Mountain and flows through the Central Range to the western plain area. The length and average slope of the main channel are 186.6 km and 31.21°, respectively. The drainage area is approximately 3157 km², making Choushui River basin the second largest drainage basin in Taiwan, with an annual runoff of approximately 6.1 billion m³ yr^− 1. Southwest monsoon and typhoons bring rainfall during the wet season from May to October. The dry season occurs from November to April, as the northeast monsoon is blocked by the Central Range. Rainfall distribution is also affected by topography, declining from mountainous areas to plains. The average annual rainfall in mountainous areas is more than 2200 mm and in plain areas is approximately 1400 mm [26].

Fig. 1

Spatial distribution of gauging stations and catchments in the Chuoshui River Basin

The Chuoshui River Basin has a unique geological environment with easily broken and weathered lithologies (Fig. 2) [27]. The downstream and midstream are divided by the mountain pass (located at the Chun-Yun Bridge, CYB, gauging station), and the downstream predominantly comprises alluvial deposits. The midstream range extends from the mountain pass to mainstream confluence with the Chenyoulan River (where the Nei-Mao-Pu, NMP, gauging is located). Most riverbanks in the midstream contain alluvial deposits, and the overall regional lithology primarily comprises sandstone, shale, and mudstone. Gravel, clay/mud, and sand are present near the mountain pass. The upstream lithology is mostly argillite, slate, and phyllite; quartzite and coaly shale appear near the midstream. Metamorphism grade of the basin lithology increases from west to east owing to the direction of orogenic movement.

Fig. 2

Spatial distribution of lithology in the Chuoshui River Basin

2.2 Data

To focus on the influence of seasonal and catchment scales on model performance and dynamic storage components estimation, we chose eight streamflow gauging stations in the Chuoshui River Basin. The geographic and hydrological information of each catchment was depicted in Tables 1 and 2. Baseflow coefficients and dynamic storage components in the dry and wet seasons were evaluated from the beginning of the wet season from 2013 to 2017. The spatial distributions of catchments and streamflow gauging stations are depicted in Fig. 1. Daily grid rainfall data were acquired from the Taiwan Climate Change Projection Information and Adaptation Knowledge Platform [28]. To estimate the water balance and effective rainfall, daily grid evapotranspiration and interception loss data were obtained from the Global Land Evaporation Amsterdam model dataset [29]. Grid data with a resolution of 0.25° × 0.25° were evaluated utilizing a series of estimation algorithms, including different land evaporation compositions [30]. Interception loss was calculated using the Gash analytical model, and actual evapotranspiration was converted from potential evapotranspiration based on observations of the microwave vegetation optical depth and root-zone soil moisture.

Table 1 Geographic information in each gauge station of the Chuoshui River Basin in Taiwan

Table 2 Hydrological information in each gauge station of the Chuoshui River Basin in Taiwan

2.3 Analytical FDC model

The analytical FDC model was developed by incorporating the rainfall–soil moisture equation proposed by Rodriguez-Iturbe et al. [31] into a model framework with rainfall-driven streamflow probability attributes [16, 17]. This model assumes that streamflow is generated from groundwater discharge and follows a non-linear recession behavior when a sequence of rainfall events increases soil moisture beyond the retention capacity [18]. Based on this assumption, Rodriguez-Iturbe et al. [31] described probabilistic dynamics of the average soil moisture utilizing a stochastic differential equation:

$$\frac{ds(t)}{dt}=-\rho \left[s(t)\right]+\xi (t)$$

(1)

where ds(t)/dt is the time derivative of the soil moisture s (mm), −ρ [s(t)] is the soil moisture loss function resulting from evapotranspiration, surface runoff, and deep percolation, and ξ(t) is the stochastic instantaneous increase caused by rainfall infiltration. Botter et al. [18] assumed that daily rainfall is a stochastic forcing of groundwater discharge generation, and applied the assumption of a non-linear recession process of daily streamflow at the catchment scale. The stochastic differential equation can be expressed as follows:

$$\frac{dQ(t)}{dt}=- aQ{(t)}^b+{\xi}^{\hbox{''}}(t)$$

(2)

where Q (mm d^− 1) is the daily streamflow, a and b are the baseflow coefficients depending on the catchment characteristics, and ξ”(t) represents a stochastic input process wherein a rainfall event provides sufficient water to generate the streamflow. The marked Poisson process with a rainfall–streamflow frequency λ (d^− 1) and the exponentially distributed rainfall depth with an average effective rainfall (which is the observed rainfall minus interception loss) on rainy days α (mm d^− 1) are assumed as the catchment rainfall input. Botter et al. [18] described the model framework as the steady-state stochastic distribution function of the daily streamflow:

$$p\left(Q,\kern0.5em t\to \infty \right)=C\left\{\frac{1}{Q^b}\exp \left[\frac{-{Q}^{2-b}}{\alpha a\left(2-b\right)}+\frac{Q^{1-b}\lambda }{a\left(1-b\right)}\right]\right\}$$

(3)

where C is a normalizing constant; a rainfall–streamflow frequency λ (d^− 1) can be obtained from the relationship between the average effective daily rainfall and average streamflow, $\overline{Q\ }$ (mm d^− 1) [32], which is calculated as follows:

$$\overline{Q}=\lambda \alpha$$

(4)

2.4 Parameter fitting

Similar to Santos et al. [19], in this study, we used maximum likelihood estimation (MLE) to calculate the baseflow coefficients a and b, providing optimal parameter fitting of the probability model for the observed data. Maximum likelihood function represents the joint probability of all observed data and can be expressed as follows:

$${\mathcal{L}}\left(b,a\right)=\prod \limits_{i=1}^Np\left(Q;\kern0.5em \textrm{b},\kern0.5em \textrm{a}\right)$$

(5)

where N is the number of data points and p(Q; b, a) is the probability density. Santos et al. [19] applied an analytical FDC model for dry seasons. Although the model performance is reduced with temporal scales from 5 to 1 year, it still simulates FDC better than REC. Additionally, all the observed data can be used in the model without data selection, which is a common process in REC.

To determine differences in the estimated baseflow coefficients for simulating the FDC, we utilized REC for baseflow coefficient estimation. REC shows that streamflow and its variation have a power-law relationship during recession period (−dQ/dt = aQ^b), characterizing the storage–discharge relationship through data fitting [7]. We selected individual recession events to conduct REC using the following criteria: (1) at least five consecutive data points for the recession event; (2) removal of data points at which the flow variation was positive or zero; (3) removal of two and one data point at beginning and the end of all recession events, respectively; (4) removal of three data points after the data point was larger than 7% exceedance probability, as defined by the FDC; and (5) removal of singular points in the data series [33, 34]. Baseflow coefficients were obtained by fitting the individual recession events using the linear least squares approach, and recession events with fitted b values < 3 were selected to ensure that the recession belonged to the aquifer drainage stages mentioned by Arumi et al. [35]. The median of the fitted b values was taken as the fixed coefficient b to re-fit a value of each recession event. The median of the fitted a was considered as the representative value [13, 36]. Additionally, we counted the number of available recession events to explore REC application in catchments and the difference between dry and wet seasons.

2.5 Performance evaluation

The probability model can compute joint probability corresponding to all the observed data, and cumulative distribution function can then be obtained by integration. To evaluate performance of the analytical FDC model at seasonal and catchment scales, we used c^KS to test the similarity between simulated and observed FDC. c^KS is an important reference value for testing whether the data distribution comes from a specific reference distribution. This is the maximum distance between the cumulative distribution functions of the simulated and the observed streamflow (F(Q) and F($\overset{\sim }{Q}$)). This value was calculated as follows:

$${c}^{KS}=\underset{x}{SUP}\left|F\left(\tilde{Q}\right)-F(Q)\right|$$

(6)

A smaller c^KS value indicates higher similarity between the cumulative distribution functions, indicating that the model has better performance. Previous studies have also applied this approach to evaluate the performance of different methods [11, 19, 22, 37]. However, it is only used to identify cumulative distribution function similarity, and its value has not been defined as a standard corresponding to model performance. Therefore, we explored applicability of the model to catchment scale and its difference in wet and dry seasons using the previous study results as a reference.

2.6 Catchment dynamic storage

Catchment dynamic storage can be regarded as groundwater storage in the unconfined aquifer, which is more easily affected by external factors. Groundwater discharge behaviors can be divided into S_d (mm) and S_i (mm). S_d is the groundwater storage that contributes to the streamflow, whereas S_i has little interaction with streamflow and changes through vertical gains or losses. Dralle et al. [25] proposed the assumption that S_T (mm) is the sum of S_d and S_i, as depicted in Eq. (7):

$${S}_T(t)={S}_d+{S}_i={\int}_0^t\left(P-Q- ET\right) d\tau$$

(7)

where P represents the daily rainfall (mm d^− 1), ET is the daily evapotranspiration (mm d^− 1), and τ is the dummy integration variable. If rainfall, evapotranspiration, and transfer between direct and indirect storage are relatively smaller than the streamflow during groundwater discharge process, the water balance will only comprise streamflow and groundwater storage changes, as depicted in Eq. (8):

$$\frac{dS_d}{dt}\approx -Q$$

(8)

Storage–discharge sensitivity function g(Q), which represents the change in direct storage corresponding to that in streamflow, can be derived as the relationship between streamflow, Q, and flow variation dQ/dt by substituting Eq. (8). S_d can be quantified by integrating the reciprocal of g(Q) using Eqs. (9) and (10).

$$g(Q)=\frac{dQ}{dS_d}=\frac{dQ/ dt}{dS_d/ dt}\approx {\left.-\frac{dQ/ dt}{Q}\right|}_{Q\gg P,\kern0.5em \textrm{ET},\kern0.5em \textrm{R}}$$

(9)

$${S}_d=\int {dS}_d={\int}_{Q(0)}^{Q(t)}\frac{dQ}{g(Q)}$$

(10)

S _i can be estimated by subtracting S_d from S_T. As it is impossible to determine boundary and condition of the catchment groundwater storage, we assumed that the initial total dynamic storage was 0. This helped us comprehend the proportion of different dynamic storage components contributing to the total amount of dynamic storage and its temporal changes under various seasonal wetness conditions.

3 Results and discussion

3.1 Baseflow coefficients

Herein, MLE was applied to fit the optimal baseflow coefficients with an analytical FDC model in dry and wet seasons, and REC was performed to estimate the coefficients to compare their results and physical information. The relationship between baseflow coefficients was also explored to determine the spatial and seasonal differences. The results are presented in Table 3.

Table 3 Baseflow coefficients with MLE and REC in dry and wet seasons

3.1.1 Estimated baseflow coefficients of two approaches

MLE results showed that the a and b ranges were 0.15–2.04 and 1.01–2.54, respectively, in the dry season. In the wet season, a ranged from 0.09 to 1.15, and b from 1.13 to 2.26. Except for BSB, LMB, and SLB (see Table 1 for abbreviations) catchments, coefficient a was lower in the dry season. Coefficient b was higher in the dry season, except in BSB and SLB catchments. Most studies [10, 38,39,40] have reported that coefficient b represents catchment recession nonlinearity, which is related to the inclination of aquifer and other hydrogeological characteristics. Bart and Hope [36] and Biswal and Kumar [41] found that coefficient a was related to antecedent catchment wetness conditions and length of the drainage system network. Coefficient a can be explicitly expressed as a function of the initial catchment storage, when coefficient b is assumed to be constant [42]. Our study did not observe a significant correlation between coefficient a and the average streamflow, which represents the catchment wetness condition. However, in SLB catchment, where the average streamflow is one order of magnitude higher than that in other catchments, there was also a similar condition for coefficient a in both wet and dry seasons. This may be largely influenced by the obvious wetness condition differences. As SLB catchment encompasses the Ming-Tan, Ming-Hu, and Sun Moon Lake reservoirs, power plant operations in this region may be related to the reservoir drawdown process. Control of hydraulic structures on rivers may decrease the frequency of low-flow events causing the model to misidentify the faster recession with a high-flow event in the SLB catchment. Brutsaert [33] also suggested that a higher coefficient a indicates the faster recession process with fewer drainage days.

REC results showed a and b value ranges of 0.001–1.82 and 1.11–2.86, respectively, in the dry season. During the wet season, the a values ranged from 0.005–0.12, and b values ranged from 1.76 to 2.36. Except for SLB, coefficient a in the dry season was higher than that in the wet season. Coefficient b was higher in CCB, YFB, NMP, and SLB in the dry season than in the wet season. Contrary to MLE results, only coefficient a demonstrated a significant seasonal difference, and the four mainstream catchments did not show a similar regional difference in the two coefficients. There was a lower average a and higher average b value in the REC than in the MLE. The catchment drainage process can be simply divided into short- and long-term states (b = 3 and b = 1.5, respectively): b = 3 represents fast drainage from the transiently saturated aquifer post initial rainfall event, and b = 1.5 represents slow drainage from the saturated aquifer post rainfall infiltration [35]. Our study found that most of the selected recession events were close to the short-term state (b = 3), resulting in a lower coefficient a. This indicates that most of the selected recession events could not completely represent discharge behavior from the aquifer. Additionally, the number of available recession events in dry and wet seasons indicated that there were less than 10 recession samples in most of the catchments over the 5-yr period analyzed (Fig. 3). Although the wet season had more recession events due to more rainfall events, the low number of available recession samples with the short-term state made it difficult to determine the average catchment discharge behavior. Faster hydrological response in catchments with small area may also be the reason for discontinuity of recession events. Therefore, REC is more suitable for long-term analysis with sufficient data.

Fig. 3

The number of selected recession events in dry and wet seasons from 2013 to 2017

3.1.2 Spatial and seasonal difference in baseflow coefficients

To explore spatial distribution of the baseflow coefficients, we plotted the relationship between the MLE baseflow coefficients (Fig. 4). The results demonstrate a negative correlation between the baseflow coefficients in both wet and dry seasons. Although coefficient a cannot be comparing owing to the difference in coefficient b, it still indicates that a depends on b. According to the coefficient characteristics mentioned above, the higher recession nonlinearity, the wetter the catchment. Among the four mainstream catchments (CCB, CYB, YFB, and BSB), there was a decrease in a and an increase in b from downstream to upstream during wet season. This demonstrates that the baseflow coefficients changed with average slope and catchment elevation, and the tributary catchments also demonstrated a similar tendency. Spatial distribution of the baseflow coefficient is uneven during the dry season. Except for the results in BSB, baseflow coefficients in the mainstream catchments had a similar tendency to the wet season and had a higher variability among the tributary catchments during the dry season.

Fig. 4

The relationship between baseflow coefficients a and b. Dry and wet seasons are indicated by red and blue symbols, respectively

According to the spatial distribution of basin lithology, hydraulic conductivity and specific yield theoretically increase from upstream to downstream, and the tributary catchments have a similar lithology, except for SLB. If the aquifer through which groundwater flows in the wet season is more extensive than that in the dry season, upstream baseflow coefficients in the wet season should reflect a lower groundwater discharge rate. This is consistent with the difference in coefficient a, and previous studies also depict that a lower a value leads to a lower recession rate [43, 44]. Therefore, seasonal differences may indicate that the groundwater discharge in dry seasons may be primarily derived from the aquifer with better drainage ability, signifying that discharge behavior also depends on the aquifer characteristics.

In most catchments, coefficient b had greater variability than coefficient a between the dry and wet seasons. Faster recession during dry season was suggested to be caused by higher evapotranspiration from the shallow unconfined aquifer [45]. However, the response of coefficient b to rainfall infiltration was inconsistent with the catchment drainage process. Coefficient b in the wet season was lower than that in the dry season, rather than closer to b = 3. According to the comparison between REC and analytical FDC model reported by Santos et al. [19], the estimated baseflow coefficient b is more similar to the model result after removing short-term stage recession events. This comparison also indicates that the analytical FDC model provides results closer to the late-state discharge process. Thus, short-term stage discharge, which is considered in REC after data selection, can be avoided.

3.2 Performance evaluation

This study utilized analytical FDC model to compute the joint probability corresponding to each streamflow data point, and the simulated FDC was obtained by integrating the probability density function. The observed and simulated FDCs for each catchment in dry and wet seasons are depicted in Figs. 5 and 6, respectively. The MLE results demonstrated that the streamflow in dry season was lower than that in the wet season, and the overall simulated FDCs were similar to the observed results, except for SLB in the dry season. This may have been caused by the sporadic low-flow probability distribution due to the higher streamflow caused by power plant operation mentioned above. This was also the reason that the lower coefficient b and higher coefficient a in SLB (Fig. 4) showed hydrological characteristics dissimilar to those of the other catchments. Simulations from the baseflow coefficients estimated by REC were evidently worse than those estimated by MLE. In REC results, low streamflow corresponded to lower probability, and high streamflow corresponded to higher probability. This was predominantly related to recession events with high streamflow owing to data selection. To comprehend the influence of baseflow coefficients on FDC, we set different coefficients as fixed values to explore changes in FDC, as depicted in Fig. 7. Although the effect of the model parameters in various climate conditions has already been discussed by Botter et al. [18], here, we primarily focused on baseflow coefficients to determine how storage–discharge relationship characterizes the simulated FDC. Higher b or lower a value leads to lower probability at a low streamflow and higher probability at a high streamflow, indicating that there are more high-flow events in the catchment. Conversely, a lower b or higher a value indicates that the catchment is dominated by a low flow. These descriptions of nonlinearity of storage–discharge relationship and streamflow magnitude are consistent with the assumptions about the short- and long-term states.

Fig. 5

Observed and simulated FDCs in dry season. (a) CCB, (b) CYB, (c) YFB, (d) BSB, (e) LMB, (f) YPB, (g) NMP, and (h) SLB catchment

Fig. 6

Observed and simulated FDCs in the wet season. (a) CCB, (b) CYB, (c) YFB, (d) BSB, (e) LMB, (f) YPB, (g) NMP, and (h) SLB catchment

Fig. 7

Schematic diagram of influence of baseflow coefficients (a) a and (b) b on simulated FDCs

To comprehend MLE and REC model performance, we used c^KS, to evaluate similarity between the observed and simulated FDCs, as depicted in Table 4. MLE results showed that c^KS range in the dry season was 0.05–0.09 and that in the wet season was 0.04–0.07. These c^KS values at the seasonal scale were expectedly higher than those reported by Santos et al. [19] at the one-year scale (c^KS ≈ 0.02). This was also close to the results of Santos et al. [22] for the summer streamflow (c^KS ≈ 0.04); therefore, the performance of MLE in our study was still acceptable. Except for YFB and NMP, c^KS values were larger in the dry season than in the wet season, indicating that the model with MLE simulated FDC superiorly in the wet season. Additionally, SLB results showed the worst performance in the dry season. This may be affected by reservoir drainage, which causes abnormally high streamflow in the dry season, and can be considered as a storage–discharge characteristic with artificial impacts. In REC, c^KS values ranges in the dry and wet seasons were 0.12–0.57 and 0.22–0.49, respectively. Except for LMB and SLB, c^KS values were lower in the dry season than in the wet season, indicating that individual recessions performed better in simulating FDC in the dry season. This is also consistent with the criterion that REC can be used only under low-flow conditions [33]. However, overall results demonstrated that c^KS values were an order of magnitude higher in REC than in MLE. Therefore, MLE results were more suitable for describing the slow and natural discharge behavior of aquifers. Nevertheless, since the model determines the probability distribution pattern, most of the higher c^KS values may have been affected by the number of data points at the seasonal scale that can be augmented to confirm whether there is a more obvious difference in the model performance between the two seasons.

Table 4 Kolmogorov–Smirnov distance (c^KS) with MLE and REC in dry and wet seasons

3.3 Dynamic storage components

After using the storage–discharge sensitivity function to estimate the S_d, we then separated S_i from S_T, which was computed using the water balance method. The dynamic storage components during dry and wet seasons are depicted in Fig. 8. The ranges of average S_d and S_i were − 11.61–7.69 mm and 0–63.89 mm in dry seasons, respectively. In wet seasons, the range of average S_d and S_i were − 16.32–24.89 mm and 0–222.83 mm, respectively. Higher S_i during the wet season indicates that more rainfall infiltration may increase aquifer storage. Most of the S_d values were negative in the dry season, representing a continuous decline in groundwater discharge. Both dynamic storage components were generally higher during the wet season, and it indicates that dynamic storage may be predominantly controlled by rainfall during the different seasons. Previous studies on groundwater level variation mechanisms have suggested that accumulated rainfall has a significant impact on groundwater level changes, and their correlation can be enhanced by setting a lower rainfall threshold [46, 47]. Chen et al. [48] also found the groundwater level fluctuation potential in the wet season was higher than that in the dry season in the Choushui River Basin through exploring the influence of dominant factors on groundwater. However, the changes in groundwater level are somewhat inconsistent with our results, and this may be caused by their estimation based on the relationship between groundwater and influencing factors rather than hydrological variables we used. It shows limitation of general water balances and importance of linking the potential influencing factors to groundwater.

Fig. 8

Average direct (S_d) and indirect storage (S_i) in (a) dry and (b) wet seasons. Dry seasons (November–April) from 2013 to 2017 and wet seasons (May–October) from 2013 to 2017 are shown from left to right

Interestingly, catchments BSB, YFB, LMB, and NMP, which had higher elevations, had the highest total dynamic storage ranges among the eight catchments. This may be related to the storage capacity of mountain aquifers. In a comparison between inclined and horizontal aquifers, Sayama et al. [2] have suggested that inclined aquifers have a larger storage capacity and smaller groundwater discharge area. Therefore, aquifer can store additional water after runoff occurs. Within the total dynamic storage, the average S_i accounted for 21 and 75% in the dry and wet seasons, respectively. These results depict that dynamic storage mostly contributes to streamflow in the dry season, and S_i accounts for more than half of the total dynamic storage in the wet season. Dralle et al. [25] found that the total amount of dynamic storage continued to reduce during the no-streamflow period in summer, which represents S_i loss due to evapotranspiration. This indicates that S_i changes are also controlled by other hydrological variables. However, whether groundwater that stores during the wet season provides baseflow for dry season is still questionable and essential for catchment ecosystem maintenance. The physical mechanism of S_i requires a more specific explanation, and there is also a high degree of uncertainty regarding water balance methods at short time scales. Hydrological models or downscale satellite gravity measurements [49, 50] can be considered to improve dynamic storage estimation in the future.

4 Conclusions

This study applied the analytical FDC model with MLE for simulating FDCs to estimate baseflow coefficients in the dry and wet seasons. Using the Kolmogorov–Smirnov distance, similarity between the simulated and observed FDCs was also determined as the model performance for exploring the difference in catchments and the two seasons and comparing them with REC. Combined with the water balance approach, the two dynamic storage components were estimated to analyze the changes and regional differences at the catchment scale. The model with MLE performed well in both dry and wet seasons and also had an order of magnitude higher performance than that with REC. Differences in the baseflow coefficients from downstream to upstream catchments in the wet seasons were consistent with the relationship between the coefficients and catchment characteristics. Seasonal difference in baseflow coefficients indicated that seasonal rainfall affected the extent of the aquifers through which groundwater flows. Additionally, a comparison with REC demonstrates that REC coefficients were mostly attributed to the short-term discharge state with faster drainage and high streamflow. Therefore, MLE results tended to characterize the storage–discharge relationship and avoid fast recession with a high-flow event. Although there was no specific trend in the dynamic storage components, catchments at higher altitudes had larger dynamic storage amounts owing to their higher storage capacity. The proportion of total dynamic storage also demonstrated an obvious seasonal difference with additional groundwater storage. Applying the analytical FDC model to estimate dynamic storage components will help improve our understanding of catchment dynamic storage response mechanisms and may offer a reference for hydrological prediction or water resource management.