Technical Note: A simple generalization of the Brutsaert and Nieber analysis

The Brutsaert and Nieber (1977) analysis is a wellknown method that can estimate soil parameters given discharge data for some aquifers. It has been used for several cases where the observed late-time behavior of the recession suggests that the water stream that is adjacent to the aquifer has nonzero depth. However, its mathematical formulation is, strictly speaking, not capable of reproducing these real-case scenarios since the early time behavior is based on a solution for which the aquifer stream has zero depth (PolubarinovaKochina, 1962). We propose a simple generalization for the Brutsaert and Nieber (1977) method that takes into consideration the depth of the adjacent water stream. The generalization is based on already available solutions by PolubarinovaKochina (1962), Chor et al. (2013) and Dias et al. (2014) and can be readily implemented with little effort. The original and proposed equations are tested against numerical simulations of the full nonlinear Boussinesq equation. A sensitivity analysis shows that the modification can have significant impact on the predicted values of both the drainable porosity and the saturated hydraulic conductivity.


Introduction
The Brutsaert and Nieber (1977) analysis (from now on referred to as BN77) has been widely used in hydrologic research to estimate aquifer parameters given some discharge data.This technique is based on "state-space"-like plots of Q × dQ/dt, where Q(t) is the aquifer discharge as a function of time.It is based on solutions for the Boussinesq equation for groundwater flow applied to a system as the one presented in Fig. 1, which shows a water channel of length L with one aquifer of length B on each side.Traditionally three solutions of the Boussinesq nonlinear partial differential equation for groundwater flow (Boussinesq, 1903) are considered for this method, which are the three solutions proposed by BN77.Solution (i) was originally given as an incomplete series by Polubarinova-Kochina (1962) (the complete set of coefficients was given by Chor et al. (2013)) for a semi-infinite aquifer that deals with early-time behavior as where the coefficients b n are given by a recurrence relation (Chor et al., 2013), D = H k 0 /n e and x is the horizontal distance.Solution (ii) is the exact solution provided by Boussinesq (1904) for a finite aquifer, which is adequate for later times.It is given by h(x, t) = F (x) where F (x) is defined from the incomplete beta function and â is a constant coefficient.Solution (iii) is the linearized solution provided by Boussinesq (1903) (also appropriate for late-time behavior) which is given by In order to obtain Eq. ( 3), the water table height h is approximated as pH in the nonlinear term of the Boussinesq partial differential equation, for linearization purposes.From the aforementioned solutions, only Solution (iii) is able to deal with nonzero water-stream depths (H 0 ) adjacent to the aquifer (of initial water table height H ). Recently, Solution (i) -from now on we call "Solution (i)" any solution for a semi-infinite aquifer where discharge is occurring -has been generalized by Chor et al. (2013) and Dias et al. (2014).The work by Dias et al. (2014) is of particular importance for us because it extends the early-time behavior to cases where the stream depth is different from zero.
Since BN77, many changes and improvements have been suggested (for detailed reviews, see Rupp andSelker, 2006, andTroch et al., 2013) but its main insight remains the same: that one should look at the rate of discharge as a function of discharge, or, mathematically (for the case of a power law), where Q is the water discharge, t is time and α and β are calibrated coefficients which can be compared to the predictions from the abovementioned analytical solutions by Polubarinova-Kochina (1962), Boussinesq (1903) and Boussinesq (1904), among many others (Rupp and Selker, 2006).
If one wishes to estimate only the soil hydraulic conductivity k 0 and the drainable porosity n e , two of the three aforementioned solutions can be used.In this work, the solutions used are the ones by Polubarinova-Kochina (1962), Dias et al. (2014), andBoussinesq (1903).However, the solution by Polubarinova-Kochina (1962) is only valid for the case H 0 = 0: it is therefore important to assess how much this assumption affects the estimate of k 0 and of n e for cases where it does not hold.
From the long list of solutions of the Boussinesq equation that are used for BN77's method, very few take H 0 into consideration (from the list of 13 equations presented by Rupp and Selker (2006), only 2 have H 0 as a parameter), so it is safe to say that the approximation of zero water level depth has not been thoroughly studied.
Although the BN77 method has been the focus of many studies for over 40 years, the subject is not, by any means, exhausted.Among recent findings is the work by Bogaart et al. (2013), which shows that, for sloping aquifers, it is possible to find a β coefficient of zero -something that until then had not been found by any other work and that was again found by Hogarth et al. (2014).Recent uses of this equation include the linking of geological and geomorphological features to hydrological behavior (Mutzner et al., 2013;Vannier et al., 2014) and the definition of good engineering practices for the robust calibration of parsimonious models (Melsen et al., 2014).
Several considerations related to the complexities of real watersheds as well as the actual physical mechanisms through which baseflow is produced and routed through the watershed raise criticism on the applicability of the BN77 recession analysis.A short, and by no means exhaustive, list of such considerations includes the effect of steep hillslopes and vertical inhomogeneity of k 0 , horizontal inhomogeneity (variation of hydraulic properties within the watershed), difficulties in the identification of α and β in Eq. ( 4) due to noisy data, geomorphological effects, etc. (Troch et al., 2013).Spatial heterogeneity appears to be particularly important and its effects were thoroughly investigated by Harman et al. (2009).This latter work shows (in a very analytical manner) that heterogeneity alone can give rise to different values of the exponent β.The values of the exponents observed in real data, however, are such that they can be explained by either the hydraulics of the aquifer or by horizontal heterogeneity (Harman et al., 2009, Fig. 9).
The usefulness of recession analysis in hydrology, however, seems indisputable, as well as the validity of the Boussinesq model in partly explaining hydrological recessions: the Boussinesq model has proved able to include realistic effects while being kept relatively simple, and remains an important tool in obtaining representative parameters for hydrological and land-surface models at the catchment scale (Pauwels and Troch, 2010;Troch et al., 2013).As such, it is reasonable to expect recession analysis and the Boussinesq model to play important roles in future progress towards improved predictive capabilities in hydrology.
It is beyond the scope of this article to explore all the considerations mentioned above.Instead, we concentrate on a single effect that has not been given much attention (H 0 = 0) and study it with a simple mathematical model that allows its importance to be assessed clearly and separately from other effects.This is in line with a systematic approach to identify inconsistencies between the theoretical models and field conditions (Pauwels and Troch, 2010, Sect. 4).Our approach using numerical solutions follows many other similar works on hydrological recessions (van de Giesen et al., 2005;Rupp and Selker, 2006;Bogaart et al., 2013) The attention in this paper is focused on the BN77 method based on Solutions (i) and (iii) (which are given in Eqs. 1 and 3) of the Boussinesq equation.We will generalize the implementation of Solution (i) with existing solutions in order to investigate the effects of the depth of the adjacent water stream into the estimation of the drainable porosity and saturated hydraulic conductivity.This generalized implementation is later shown to considerably improve the estimation of the hydraulic conductivity and drainable porosity in numerically generated data.This result suggests that this new formulation of the BN77 analysis could potentially be useful for obtaining k 0 and n e in man-made drainage systems and improving simulations of drainage and water table dynamics from hypothetical hillslopes and for better understanding single hillslope processes, where the BN77 analysis is more likely to succeed (Troch et al., 2013, Sect. 4.3).

Generalization of the early-time equation
Let ξ denote the Boltzmann variable for the one-dimensional Boussinesq equation (Chor et al., 2013;Dias et al., 2014), where D = H k 0 /n e , and φ denotes a normalized water table height, where h(x, t) is the water table height, x is the horizontal distance from the water stream and t is time.Under the above change of variables, the Boussinesq equation is reduced to the dimensionless ordinary differential equation together with the boundary conditions φ(0) = φ 0 and φ(∞) = 1.Due to the second boundary condition, the solution is only valid for the initial phase of aquifer drawdown.
For φ 0 = 0, as already noted, the solution by Polubarinova-Kochina (1962) suffices for the BN77 analysis; for φ 0 = 0, a series solution of the form has been proposed by Dias et al. (2014), with a recursion relation for the a n 's.An important result in that work is an empirical equation, fitted to numerically obtained values of a 1 in the series above, for the value of ψ 0 , defined below.This is given as Eq. ( 16) in the present work.
Let us also define which we apply to Darcy's law, along with Eqs. ( 5) and ( 6) to obtain where q(x, t) is the flow rate per unit width at any point x of the aquifer.Since we are interested in the aquifer-stream interaction, we set x = 0, which produces where ψ 0 ≡ ψ(ξ = 0) and φ 0 ≡ H 0 /H .The value of ψ 0 , as far as we know, cannot be obtained analytically and is generally obtained numerically or by means of approximations: its calculation will be dealt with later.For now, it suffices to note that ψ 0 is a function of φ 0 as given above.
Writing dQ/dt = −α 1 Q β 1 , where the subscript 1 indicates the early-time solution, and Q = 2L q is the flow per unit length taken over the total length (L) of the tributary and main channel sections upstream from the gaging station, with q as in Eq. ( 11), yields β 1 = 3 and Equation ( 12) is generally used with the assumption of H 0 = 0, which yields ψ 0 (0) = 0 ≈ 0.6642, which (substituting back into Eq.12) gives the well-known Eq. (18b) of BN77.
However, often the value of H 0 is not small enough in comparison with H in order for this approximation to be valid (Munster et al., 1996;Serrano and Workman, 1998;Barlow et al., 2000;Peterson and Connelly, 2001;Langhoff et al., 2006;Ha et al., 2008;Sena and de Melo, 2012).In these cases the misplaced assumption could lead to biased estimates of k 0 and n e .These latter errors depend not only on the determination of α 1 but also on the late-time equations chosen and on the determination of the constants for that solution.
Evidence that the water depth of the adjoining stream is not negligible can be found (for example) in the work by Brutsaert and Lopez (1998), where the late-time data showed a decay with β 2 ≈ 1, which in fact indicates that the watershed analyzed has a ratio H 0 /H close to one (we use the subscript 2 to indicate the late-time solution).Indeed, for φ 0 = 0, the exact analytical solution provided by Boussinesq (1904), which is valid for late times, gives β 2 = 3/2 (Brutsaert and Nieber, 1977), whereas numerical solutions of the Boussinesq equation (Kan, 2005) show that β 2 varies from 3/2 down to 1 as H 0 /H varies from 0 to 1.
Comparison between the saturated hydraulic conductivity k 0 as estimated by the original BN77 analysis (dotted line and symbols) and the method proposed here (solid line) -both normalized by the real value k 0 used to numerically generate the discharge data and varying against the normalized water height at the origin φ 0 .

Comparison between both approaches
We dedicate this section to the estimation of the errors that arise by assuming that the stream depth H 0 is zero.For that purpose we take as a late-time equation the solution of the linearized Boussinesq equation, given by Eq. ( 3), which predicts β 2 = 1 and where A is the area of the watershed, approximated by 2B L. Solution of Eqs. ( 13) and ( 12) gives, for n e and k 0 , and In this formulation we assume both ψ 0 and p to be functions of φ 0 = H 0 /H , so we have ψ 0 (φ 0 ) and p(φ 0 ), as was previously emphasized.We also assume that p(φ 0 ) = (1 − p 0 ) φ 0 + p 0 , where p 0 = 0.3465, based on the fact that p = 0.3465 for H 0 = 0 (Brutsaert and Lopez, 1998).Setting H 0 = 0 (and therefore φ 0 = 0) in this model will yield exactly the same equations as presented by Brutsaert and Lopez (1998).
To obtain ψ 0 (φ 0 ) we use the approximation provided by Eq. ( 14) of Dias et al. (2014), since it is sufficiently accurate and simple to program, viz.for the a n 's in Eq. ( 8) has been obtained, the values of the a n 's still cannot be obtained analytically, essentially because the series' radius of convergence is limited so that the boundary condition φ(∞) = 1 cannot be imposed analytically.Instead, they must be obtained numerically with the aid of numerical solutions of Eq. ( 7).The coefficients above have been obtained in Dias et al. (2014) by curve fitting with a large number of numerical solutions.
In order to compare both approaches, we have solved numerically the original Boussinesq equation (in x and t) to model the system depicted in Fig. 1 and generated synthetic discharge data for different values of H 0 < H .We then applied the original BN77 method to these data, as well as the generalized method we propose here.With this analysis we can quantify the error of both methods in order to determine their accuracy.
Figures 2 and 3 show the results for k 0 and n e , respectively, for increasing values of φ 0 , plotted against the true values k 0 and n e .As can be seen, the k 0 estimate using the original equations remains close to the true value up to φ 0 = 0.4 approximately.Furthermore, there is a more or less linear trend in n e estimated with the original equations all the way from φ 0 = 0.Both differ considerably from the truth for large values of φ 0 (φ 0 > 0.4 for k 0 and φ 0 > 0.2 for n e , approximately).On the other hand, our Eqs.( 15) and ( 14) give estimated values of k 0 and n e that differ very little from the true ones for the whole φ 0 range, and as such represent a considerable improvement over the original equations.
The small kinks between φ 0 = 0.7 and φ 0 = 0.8 are an artifact of the choice of the range of the streamflow Q for fitting α 1 and α 2 used in the recession analysis.This (to the best of our knowledge) is still a subjective part of the BN77 analysis: the ranges were chosen to fit the recession plots dQ/dt × Q reasonably well, but they were not "fudged" to "optimize", in any way, the estimated k 0 and n e .

Conclusions
We have given an expression for early time aquifer discharge that generalizes the broadly used Eq. ( 18) of Brutsaert and Nieber (1977) for cases where H 0 is not small enough compared with H to make φ 0 = 0 a valid approximation and compared the results to the original BN77 method.The main motivation for this approach was to investigate the effects of this assumption on the determination of the saturated hydraulic conductivity k 0 and drainable porosity n e .This generalization, given mainly by Eq. ( 12), is easily applicable and requires virtually no change in the original theory presented by BN77.
The comparisons presented in Figs. 2 and 3 suggest that the BN77 estimates of the hydraulic conductivity k 0 underestimate the "true" numerical value more and more as φ 0 increases.This under-estimation is particularly large when H 0 is above 40 % the value of H (φ 0 > 0.4).The same behavior is observed with the estimation of the drainable porosity n e , with the under-estimation being specially large for H 0 greater than 20 % the value of H (φ 0 > 0.2).On the other hand, the method presented here deviates very little from the "true" numerical values in both cases and for any value of φ 0 .We consider the errors using the original BN77 analysis for both cases to be large enough that the water stream depth should be considered as a variable when using BN77 to estimate these parameters.

Figure 1 .
Figure 1.Schematic of a watershed of simple geometry during a hydrologic recession.

Figure 3 .
Figure3.Comparison between the porosity n e as estimated by the original BN77 analysis (dotted line and symbols) and the method proposed here (solid line) -both normalized by the real value n e used to numerically generate the discharge data and varying against the normalized water height at the origin φ 0 .