New Analysis Method for Continuous Base-Flow and Availability of Water Resources Based on Parallel Linear Reservoir Models

: Water ﬂows in the hydrosphere through a tangled and tortuous labyrinth of ways that is the hydrological cycle. Flow separation models are an attempt to group such complexity of paths into a few components of ﬂow and storage so as to reﬂect the overall behaviour of a basin. A new method of analysis and separation of ﬂow components, based on equations of dynamic relations between Linear Reservoirs connected in Parallel (PLR models), is developed in this article. A synthesis of models based on mathematical ﬁlter equations is carried out in order to make comparisons with the proposed model. Reference is also made to the methodology of adjustment and calibration of the PLR models based on the recession curves of the real hydrographs. The models are tested with the continuous register of a basin located in the northeast of Spain. The simulations are carried out with two reservoir models (2R models), three reservoirs (3R models) and with a mathematical ﬁlter model to compare the results. With the results of the models, ﬂow duration curves (FDCs) and storage duration curves (SDCs) were elaborated, thus allowing assessment of the origin of the water resources of the basin, a guarantee of their regulation and availability, the dynamic storage in the catchment, residence times and other features.


Introduction
The generation of runoff is surely the most visible process of those that make up the hydrological cycle. Discharge is an integrated signal over a larger region, unlike precipitation, which is a point measurement that suffers from spatial sampling error and gauge undercatch [1]. In rivers with big watershed areas, rivers with a low fluctuation regime throughout the year and for estimations of the continuous or annual contribution of a river, better efficiencies are obtained with the direct measurement of the flow using gauging stations; however, this is not always true, especially in rivers that show a very irregular regime [2]. The processes that make up the hydrological cycle maintain marked connectivity relations. There are many processes that depend on the generation of runoff, such as erosion and sedimentation [3][4][5][6]. The runoff itself circulates through different interconnected environments, such as the River-Aquifer-Spring relationships [7]. separation is established, quick and slow; but with PLR models there may be intermediate components that sometimes represent the actual reservoirs of the basins well; (3) Unlike other methodologies, with this one, if a flow component or the total discharge is known, the DRE allows the deduction of the other flow components as a single solution, without having to make estimates of flow sharing. This last feature makes PLR models useful for the simulation of certain hydrological processes, such as rainfall-runoff transformation and routing, at different parts of the basin. The separation of the flow components-the aim of the present article-will be framed within this last process.
The formulation originates from two classic equations-the discharge or storage equation and the continuity or water balance equation. In this article, an adaptation of this formulation is made to be applied to the separation of the base-flow from the continuous time series of total streamflow.
The main objective of this work is to present the PLR model as suitable for carrying out the base-flow separation from the time series of total discharge. For this purpose, the methodological development for the PLR modelling is presented based on the recession curves of the real hydrographs. Subsequently, the application of the model to a basin in the northeast of Spain is presented. In this work, the results of the PLR models are compared with the results of other models already diffused and accepted as those based on the Mathematical Digital Filtering Method. The work was performed using software developed by the authors of this article in the Department of Earth Sciences at the University of Zaragoza. The SHEE (Simulation of Hydrological Extreme Events) package is an adaptation of traditional hydrological models to DEMs and databases. It has resulted in several publications, including those related to hydrology [2,[39][40][41][42][43][44][45][46], and it uses powerful libraries (e.g., OpenGL, GDI, GDAL, Proj4) for the management and display of DEM and datasets. Additionally, its interface provides rapid and high quality OPENGL graphics, in both RASTER and VECTOR formats.

Mathematical Digital Filtering Method
Mathematical Digital Filtering is a recursive digital filter method commonly used in signal analysis. Reference [47] proposes Equation (1), subject to F t ≤ Q t . Reference [48] adapted the method for stream flow partitioning, making the assumption that low frequency base-flow could be distinguished from high frequency high flows and proposed Equation (2). References [49,50], with the collaboration of other authors, provide a fruitful trajectory on algorithms based on Digital Filters, with Equations (3) and (4) being two of the most widespread. The base-flow index (BFI) is the ratio of the base-flow volume (calculated from F t ) to the volume of streamflow (Q t ) over a specified period. Reference [51] proposed the general form (Equation (5)) of a digital filter considering a digital filter parameter and BFImax (maximum value of long term BFI over a specified period).
where F t is the quick-flow response at the tth sampling instant, Q t is the original streamflow at the tth sampling instant and α is the parameter that enables the shape of the separation to be modified.
where B t is the base-flow at the tth sampling instant. The filter is generally run through multiple times per dataset. For example, three passes are commonly used-forward, backward and forward.

Parallel Linear Reservoir Models
In the PLR models, for a given instant, the dynamic relations between the different deposits that compose a model are controlled by some of the forms of Equation (7) [39], where nr is the number of reservoirs (at other times they are known as cells and as deposits) that make up the model (usually works with models of 2 or 3 reservoirs); Q i is the discharge of the reservoir I, and S i is the storage. For a model of two reservoirs, one is the quick-flow and the other is the base-flow.
Q 0i and α i are parameters of the model that are determined from the recession curves of the hydrographs. The α i parameter is called a depletion or recession coefficient and is also known as the response factor. Q 0i is the flow for a given instant, being the same instant for all the Q 0i . To be precise, there are nr-1 equations and since we are trying to find nr unknowns, another equation is needed-the mass conservation equation (Equation (8)).
Once the model is calibrated, that is, once the parameters Q 0i and α i are already known based on the streamflow records (Q t ) of the time series, figuring out the Q i values of each reservoir through Equations (7) and (8) is needed, thus separating the quick-flow and the base-flow. Therefore, it is possible to propose a unique equation dependent on the discharge of one of the reservoirs (Equation (9)) with a single unknown, Q 1 .
Since it is a nonlinear equation, it can be solved by successive iterations by applying Taylor's development to linearize the equation and solve it using the Newton-Raphson method (see Appendix A). Once Q 1 is obtained through the set of Equation (7), the remaining Q i values can be obtained.
For a given instant, once the distribution of the discharge of each reservoir is known, Equation (10) can determine the storage (S) or natural water reserves of the entire basin and the available water reserves in each reservoir.

Calibration of the Parallel Linear Reservoir Models
In order to be able to use the methodology described in the previous section, it is necessary to establish the Linear Reservoir Model, that is, to adjust the parameters Q 0i and α i to the actual recession curves of the time series of total streamflow.
The recession curves are therefore hydrograph sections conditioned by the release of water from natural reservoirs, that is, they are the result of the sum of the multiple contributions from natural reservoirs. To calibrate the model (i.e., to obtain the Q 0i and α i parameters), segments from the recession curve are selected from an actual hydrograph. These segments can be individually or collectively analysed to understand the processes that generate the flow. To obtain a solution for these segments, a graphical approach is traditionally adopted. However, recently analysis has involved the definition of an analytical solution or a mathematical model that can properly fit the recession segments. Reference [39] developed a completely automatic method that is based on adjustment via the least squares method, combining the Reference [12] equation with Equation (7) of this paper.

Case Study of the Continuous Streamflow Time Series of a Catchment
The application of the described methodology was carried out with the streamflow time series (  The main materials that constitute aquifers in the Alcanadre river basin are the limestones and conglomerates of the Miocene and Eocene that appear in the Sierra de Guara and the conglomerates of the marginal zone of the Ebro basin While being less important than the previous aquifers, it is possible to emphasize the existence of quaternary aquifers created from alluvial deposits of the rivers and the extensive and low-thickness deposits of gravels and clays that constitute glacis.
In Figure 2, the catchment is presented with its water divide and a geological synthesis, the altimetry distribution and the isohyets of annual precipitation. As can be seen, the northern zone is the main zone for the recharge of aquifers, due to it being the one with the highest precipitation and where calcareous formations-which represent the origin of the main aquifers of the basin-outcrop. In Table 1, some characteristics of the basin are presented. In relation to this basin, there are important antecedents among which we highlight the work of Reference [52], which performed an analysis of the recession curves of the rivers Alcanadre and Guatizalema to perform a hydrodynamic characterization of carbonated aquifers. This is related to the dynamic volume of the basin with the streamflow of the rivers, which allows us to make its corresponding streamflow recession theoretical curve.

Setting PLR Models
To set the Reservoir Models, the 14 recession curves shown in Figure 3 were selected. The adjustment methodology applied resulted in the parameters shown in Tables 2 and 3 for the 14 curves. The adjustment has been made for two models, a two reservoirs model (2R model) and a three reservoirs model (3R model).    Figure 4 shows the parameters obtained in the 14 recession curves, expressed in logarithmic scales, −Ln(α) versus −Ln(Q 0 ). In the case of the 3R model, the parameters of the 14 recession curves are aligned. Each model is the best fit for a particular recession curve but it is necessary to establish a generic model applicable to the complete time series. The mean of all the parameters of each model is the result yielding the best fit by the least squares method and approaching the regression line. Tables 2 and 3 show the parameters for the generic 2R model and the 3R model, respectively. The squared Pearson's correlation coefficient (R 2 ), with respect to the linear regression line, is 0.789 for the 2R Model and 0.934 for the 3R model.

Flow Separation with PLR Models
With the Q 0 y α parameters of the generic models described in Tables 2 and 3, flow separation in each reservoir for the continuous time series of the stream gauge A033 was performed. The series extends temporarily from 1987 to 2013, with data recorded every 15 min. In Figures 5 and 6, six stretches of the hydrograph of the time series are represented.  The base-flow predicted by the 3R model is smaller than that of the 2R model but the 3R model takes into account an intermediate reservoir that has a notable influence on the regulation of the basin. From a purely mathematical point of view, both the 2R and 3R models can be adjusted to the real recession curves with great performance, but in terms of conventional concepts of base flow separation, the 3R model appears to be more plausible than the 2R model in this particular catchment due to its characteristics and especially due to the development of a karst area in the head zone. This basin could be described with three flow components: a fast one, corresponding to the superficial flow, which possesses great relevance during the river floods; an intermediate flow, which corresponds mainly to the circulation of the larger karstic galleries; and a third, slow flow, which would correspond to the circulation through the joints, fissures, smaller diaclases and to the aquifers developed in granular soils in the alluvial terraces of the river.

Flow Separation with Mathematical Digital Filter
Simulations were performed with the Mathematical Digital Filter models of Equations (1)-(5) and it has been verified that the differences between them are not as relevant as the differences between the PLR models. This is the reason why only one single Digital filter model was represented, the one corresponding to the equation in Reference [49], due to it being one of the most widespread. Figure 7 shows the flow separation results for the six stretches. The digital filter model only allows separation of the base stream without considering an intermediate flow.

Efficiency Analysis
In order to compare the proposed reservoir models, a benchmarking analysis of the base flow obtained with the classic filter model and with the 2R model has been carried out. The results of this study for the six stretches are in Table 4.
Several indices have been used: the Nash-Sutcliffe efficiency coefficient [53], the index of agreement (IOA), the correlation coefficient of the Pearson product between the observed and calculated data (R), the coefficient of determination (R 2 ), the root mean square error (RMS, RMSE or EMC), the ratio of absolute error of peak flow (RAE), the Nash-Sutcliffe efficiency coefficient with logarithmic values (lnNSE), the Relative efficiency criteria of the Nash-Sutcliffe coefficient (NSErel), the Relative efficiency criteria of the index of agreement (IOArel), the weighted coefficient of determination (wR 2 ) and the percent bias (pBIAS). The meaning of these parameters and their equations can be found in Reference [39] and for pBIAS in Reference [54].
On the other hand, in order to compare these results with similar studies, Reference [55] made a benchmarking with the NSE and pBIAS indices, with the results obtained with baseflow separation methods, versus physically based models. The study is carried out with eight main scenarios where two kinds of terrain, sandy and silt-sandy, three different states of the groundwater level and the presence of several water pumping scenarios from the aquifer are contemplated. The results obtained have a high range of variation. For NSE, they obtained values between 0.1 and 0.9. In our study, we obtained 0.344 for stretch 1 and values of around 0.8 for the other stretches. With pBIAS, these authors obtained values between −16% and +34%. In our case, we obtained negative values in stretches 3, 4 and 6, with a minimum value of −5.1% and a maximum value of 25.5% in Section 1.

Discussion
The equation of Dynamic Relations between deposits (Equation (7), [39]) allows for the establishment of continuous streamflow separation models in two or more flow components, quick-flow, base-flow and intermediate flows for other cases. These models are calibrated with very few parameters (Q 0i and α i ), thus they are parsimonious models, which fit the recession curves of the real hydrographs of the streamflow time series. Figures 5 and 6 show the flow separation results for the 2R and 3R models and Figure 7 describes the results for the filter model of Reference [49].

Flow Duration Curves
Based on the complete time series and the flow separation results, Figure 8 shows the Flow Duration Curves (FDCs) for the three models, which are cumulative frequency curves that show the percentages of time during which specified discharges were equal or were exceeded in given periods (exceedance probability). FDCs represent the flow characteristics of a stream throughout its range of variation [38]. FDCs can also be determined for the separated stream components, thus allowing the assessment of the characteristics of each flow component and therefore the characteristics of the natural deposits of the basin assigned to each stream. In order to establish the FDCs, the streamflow values were ordered from highest to lowest and accompanied by the flow components associated with each streamflow value. As can be seen, for the mathematical filter model, the values of the associated flows do not keep a good correlative distribution in terms of order.
On the other hand, flow duration curves for streamflow, base-flow, quick-flow and intermediate flow all follow a decreasing distribution. The differences between the curves of each flow component are evident. The quick-flow duration curve (QFDC) decreases very quickly compared to the base-flow duration curve (BFDC), which decreases more slowly. The QFDC has very high values at the beginning but most of the time presents values lower than those of the BFDC. The curves of each component are cut at certain points-at one point for the 2R model and at three points for the 3R model. Figure 8D shows the cut-off points for the 3R model in detail and Table 5 presents the cut-off values for the two PLR models. Within a 0.13% to 3.94% range (0.46 to 14.39 days/year), the intermediate reservoir is the one contributing the most to the total streamflow of the river. Only at 0.46 days/year and coinciding with flood episodes of the river, it is the quick reservoir that contributes the bigger discharge. However, it is the base-flow reservoir that contributes more for most of the time-about 350 days/year.

Contribution of Each Flow Component to the Total Streamflow
If FDCs are integrated, the water volume supplied by each flow component to the total streamflow can be obtained. Table 6 shows the results expressed as average annual volume and as a percentage of each reservoir over the total streamflow. For the 3R model it follows that, for the Alcanadre river at the A033 gauging station, the annual average streamflow is 117 hm 3 /year, of which 8 hm 3 circulate through the quick reservoir, 33 hm 3 circulate through the intermediate reservoir and 76 hm 3 flow through the base flow reservoir. With all these data, the importance of the slow reservoir for two aspects is highlighted: the first aspect is regarding the annual volume of water circulating for the basin-76 hm 3 -which represents 65% of the total streamflow, while the second aspect refers to the capacity of regulation of the water resource available in the basin and distributed during most of the time, which for 350 days/year represents the main contribution to surface water.

Availability of Water Resources in the Catchment and Residence Time
In addition to the FDCs discussed above, Equation (10) also provides Storage Duration Curves (SDCs), which represent cumulative frequency curves that show the percentage of time during which specified availability of water resources were equalled or exceeded in given periods (exceedance probability). Figure 9 shows these curves for the two PLR models and, from them, interesting information can be extracted. The storage in the slow reservoir has two well-differentiated phases, one characterized by the other reservoirs being active and the other one, a longer one, when the other reservoirs have no influence.  Table 7 shows the maximum, minimum and mean values of discharge and storage for each reservoir. The volume of water stored in a watershed that is susceptible to influence the flow of the river is known as dynamic storage (DS) and can be estimated as the difference between the maximum and minimum storage. For the 3R model, the DS is 23 hm 3 , with a distribution of 7.1, 7.3 and 8.7 hm 3 respectively for quick, middle and slow reservoirs, respectively. Water availability in each deposit varies with time. During flood episodes, reservoir recharge occurs and when the rains cease, the reservoirs discharge at different velocities. The total storage for each time step is different for each model because it is the sum of the storage of each reservoir and is calculated by Equation (10) from the discharge (Q i ) of each reservoir.
The mean storage (MS) represents the average availability of the water resource along the complete time series. The 3R model has a mean storage of 0.01, 0.13 and 2.91 hm 3 for reservoirs 1, 2 and 3, respectively, which again indicates that the main water guarantee comes from the slow reservoir.
The mean discharge (MD) represents the average flow throughout the year and is calculated as the continuous flow for a year required for discharging the annual volume. Note that mean storage and mean discharge do not correlate between models, with the maximum value of the mean storage for the slow reservoir of the 3R model being 2.91 hm 3 So, for the 3R model, the residence time is 0.24, 1.45 and 14.03 days for the reservoirs quick to slow, respectively.
The storage for each reservoir can also be calculated from the FDCs, as shown in Table 8. Note that at the cut-off points, in two of the tanks, the same flow rate is discharged but the storage differs considerably.

Conclusions
In this work, a new method for base-flow separation based on Parallel Linear Reservoir models is proposed. The Dynamic Relations equation sets up a single and optimal solution for the parameters that govern the Parallel Linear Reservoir models. With these models, more than two flow components can be established to assess the water routing through the basin. The Flow Duration Curves allow evaluation of the dynamic relations between flow components that circulate through the labyrinths of the hydrological cycle of the basin. The Storage Duration Curves allow assessment of the temporary availability of the water resources in the basin.
This article demonstrates the suitability of the PLR models and the equations that describe them, to elaborate models of continuous development for the separation of different flows from the total streamflow recorded in gauged stations. Based on the assessment of the results of the PLR models, through flow duration curves (FDCs) and other techniques, evaluations can be made about the dynamic relations between different components of the flow, the recharge of aquifers and the state of the water resources of the basins.
The developed methodology has interesting specific characteristics and a series of properties that makes it suitable for such purposes, among which we can highlight the following: − Calibration of model parameters, based on actual recession curves recorded in the hydrographs of gauging stations. The proposed formulation allows adjustment of the parameters of the models with a single optimal solution, derived from a least squares adjustment, unlike other techniques in which it is necessary to manually intervene in the adjustment. This makes the results of the PLR models conform to the real hydrographs, especially in the recession stretches, as demonstrated in the study of actual episodes carried out in Reference [39]. − Models with more than two flow components can be established. In the example proposed for the Alcanadre River in Spain, three components better reflect the actual response of the basin, in the upper zone of which a large karst system is developed. In this area, the main recharge of fracturing aquifers developed in the carbonate formations (ages between the upper Cretaceous and the Eocene) that outcrop in the more elevated parts of the catchment takes place. − In addition to the separation between the different flow components, the PLR models with the proposed formulation allow the investigation and evaluation of the dynamic relations between the different components and also the drawing of conclusions about their state in each moment, the volumes stored in each reservoir, the discharges over time, the water reserves in each reservoir and the entire basin, the detractions of the aquifer flow and its recharge, the average residence time of the water volume, the volume of contribution of each type of reservoir, and so forth.
Special mention must be made of the state of water resources in the basin. In the example of the Alcanadre River in Spain, relevant information has been obtained, which can be used for watersheds in other places by following the methodology set out in this paper. For example, the importance of the slow reservoir in terms of water supply, which accounts for 65% of the total contribution, is considered quantitatively and, in terms of water reserves and water-regulation guarantee, is equivalent to a reservoir with 26 hm 3 of capacity.
were p i = α i/α 1 (A2) and m i = Q 0i Q 01 p i (A3) The first derivative of the function is: And, therefore, for iteration k, the corrected value is as follows: and so, through successive iterations, we arrive at the solution.