A time-series analysis framework for the flood-wave method to estimate groundwater model parameters

The flood-wave method is implemented within the framework of time-series analysis to estimate aquifer parameters for use in a groundwater model. The resulting extended flood-wave method is applicable to situations where groundwater fluctuations are affected significantly by time-varying precipitation and evaporation. Response functions for time-series analysis are generated with an analytic groundwater model describing stream–aquifer interaction. Analytical response functions play the same role as the well function in a pumping test, which is to translate observed head variations into groundwater model parameters by means of a parsimonious model equation. An important difference as compared to the traditional flood-wave method and pumping tests is that aquifer parameters are inferred from the combined effects of precipitation, evaporation, and stream stage fluctuations. Naturally occurring fluctuations are separated in contributions from different stresses. The proposed method is illustrated with data collected near a lowland river in the Netherlands. Special emphasis is put on the interpretation of the streambed resistance. The resistance of the streambed is the result of stream-line contraction instead of a semi-pervious streambed, which is concluded through comparison with the head loss calculated with an analytical two-dimensional cross-section model.


Introduction
The development of methods to estimate aquifer parameters from stream-aquifer interaction dates back to the 1960s and the early application of computers in hydrology (Cooper and Rorabaugh 1963;Pinder et al. 1969;Venetis 1970). The approach proposed at that time, referred to as the flood-wave method, is similar to a pumping test, as the groundwater head in an aquifer is perturbed by a single stress, in this case a flood wave in a stream adjacent to the aquifer. The aquifer diffusivity is obtained by fitting a simple equation for stream-aquifer interaction to the observed heads. This equation fulfills the same function as the well functions of pumping tests. Hall and Moench (1972) refined the method by using convolution integrals to relate stream stage fluctuations and head fluctuations. Later, Moench and Barlow (2000) extended the method by adding equations for a set of different stream-aquifer configurations. Alternatively, groundwater head response to a time series of stream stage fluctuations can be obtained analytically by replacing the time series of observed stream stage by a series of basis splines (Knight and Rassam 2007;Rassam et al. 2008).
A limitation of the flood-wave method is that it is applicable only to situations where head fluctuations can be clearly related to river stage fluctuations (Ha et al. 2007).
Electronic supplementary material The online version of this article (doi:10.1007/s10040-016-1436-5) contains supplementary material, which is available to authorized users. In many cases, however, this is not possible as fluctuations due to other stresses, like recharge and evaporation, interfere with fluctuations due to stream stages variations. To solve this issue, the influence of each stress needs to be identified separately. This is where time series analysis can improve the flood-wave method.
The objective of this paper is to embed the flood-wave method into a time-series-analysis framework in order to derive aquifer parameters for use in distributed groundwater models. The framework is the method of predefined response functions (Von Asmuth et al. 2008), in which a specific response function (also referred to as a transfer function) is chosen for each stress. Each function is able to simulate the head response due to an impulse of a specific stress. Convolution of each response function with the corresponding stress time series results in the separate fluctuations caused by each stress, where it is assumed that the system's response is linear. The method of predefined response functions has recently been extended to simulate nonlinear reactions of the phreatic water table in Australia (Peterson and Western 2014;Shapoori et al. 2015a, b, c). An evaluation of the method using synthetic data was presented by Shapoori et al. (2015a, b, c). Another extension of the method concerns the estimation of aquifer parameters from time series analysis in the vicinity of well fields (Obergfell et al. 2013;Shapoori et al. 2015a, b, c).
Typically, the selected response functions do not depend on physical parameters. For example, a scaled gamma distribution function is commonly used as the impulse response function for groundwater recharge. The novelty of this paper is two-fold-first, an analytical groundwater model is used as the predefined response function similar to the functions used in the flood-wave method; second, the flood-wave method is placed in the framework of time series analysis. The resulting approach is an extension of the flood-wave method in the sense that it is applicable to situations in which other time-varying stresses than stream stage variations have a significant effect on head fluctuations. This paper is organized as follows. First, the method of time series analysis by predefined response functions is reviewed and it is explained how the flood-wave method can be placed in a time series framework. Next, a description of the hydrogeological situation of the field site is given for which response functions are developed. The time series model is fitted to data collected near the Dutch lowland river 'Aa', and aquifer parameters are estimated. These parameters are then entered into a numerical distributed groundwater model to evaluate their adequacy as parameters estimates. The physical significance of the parameter values is discussed, with a special emphasis on the interpretation of the resistance of the streambed.
Review of time-series analysis with predefined response functions

Response functions
In this paper, the flood-wave method is placed in a time-seriesanalysis framework. Time series analysis is performed with the method of predefined response functions (Von Asmuth et al. 2002). Transfer functions, a term widely used in system theory and time series analysis, can be considered as synonymous to response functions. Similar to linear systems theory (Hespanha 2009), output signals are obtained by convolution of response functions with input signals. Response functions are mathematical expressions relating input and output signals (Box and Jenkins 1976). In this paper, groundwater systems are approximated as linear in the sense that output signals are proportional to input signals. Hydraulic stresses like precipitation, evaporation, river stage variations, and pumping are the input signals and head fluctuations form the output signal. Conditions for when the approximation of linearity is reasonable are reviewed in Barlow et al. (2000).
A time series of head fluctuations φ(t) at a specific point in space can be obtained by convolving a stress time series p(t) with the corresponding impulse response function θ(t): where t is time. In this paper, φ(t) is used for head fluctuations caused by one specific stress, while h(t) is used to refer to the head fluctuations caused by the superposition of all stresses. The dimension of θ(t) is determined by the dimension of the stress so that the product p(τ)θ(t−τ)dτ has the dimension length, like heads-in contrast to linear system theory, where transfer functions are dimensionless (Hespanha 2009). Note that the dependence of the response function on spatial coordinates is omitted in this notation. The response function can also be interpreted as the weighting function in a moving average process (Olsthoorn 2008). As a comparison, in runoff hydrology, the familiar unit hydrograph is the response function relating precipitation (the input signal) to stream discharge (the output signal).
The response function of precipitation represents the passage through the unsaturated zone, followed by a recession curve describing the subsurface drainage of the infiltrated water (e.g., Besbes and de Marsily 1984). A first approximation for the response function of evaporation is the response function of precipitation multiplied by a negative scale factor. Alternatively, evaporation can be attributed its own response function describing; for example, how the root zone reacts to a drought period (Peterson and Western 2014). The response functions for river stage variations and pumping represent the propagation of the head change from the river or the pumping well to a point in the aquifer.

Discrete inputs and continuous transfer functions
In this section, it is described how time series of groundwater heads are modeled given discrete time series of stresses and continuous transfer functions. The unit step function s(t) is obtained from the impulse response function θ(t) as The step function has the dimension of length per dimension of stress.
The function is called the block response and represents the response to a unit stress applied from t = 0 to t = Δt. In this paper, the block function is used as the response function of a given stress. Time is discretized in stress periods, where Δt i is the length of stress period i. Stress p i is constant over stress period i from t = t i to t = t i + Δt i . Since the system is approximated as linear, the head at time t j can be obtained by summing the effects at time t j of all past stress periods: where For application of the extended flood-wave method in this paper, the heads h(t) in the aquifer are obtained as the following sum: where h(t) is the head, d is the drainage base which is defined as the head that is reached when all stresses are zero, and φ p (t), φ e (t), and φ s (t) represent the contributions of precipitation, evaporation, and stream stage respectively. n(t) represents the residual time series defined as the difference between observed and simulated heads. If the characteristics of the residual time series substantially depart from white noise, modeling the residual is recommended (Von Asmuth and Bierkens 2005). In this paper, an exponentially decreasing noise model is adopted.

Field site
The field site is located in the area managed by the Dutch Water Board Aa and Maas in the southeastern part of the Netherlands (Fig. 1). Piezometers were placed by the Water Board, perpendicular to the River Aa, as part of a larger monitoring program of groundwater levels. The Aa is a small, 67km-long lowland river, with an average flow of 11 m 3 s −1 at its mouth.
The field site is situated near the eastern edge of the Dutch Central Graben. The edge of the graben is a fault zone of low permeability, referred to as the Peel border fault zone. The graben is subsiding since the beginning of the Oligocene (ca 25 million years ago) and is filled with sediments over a thickness of approximately 2,000 m. Regional bore logs from the Dutch Geological Survey in the vicinity of the field site suggest that a clay layer is present at a depth of approximately 30 m bgl. This clay layer belongs to the fluvial formation of Waalre, deposited by the Rhine about 2 million years ago. The clay layer is approximately 1 m thick and can be considered as the impermeable base of the hydrogeological system.
The system above the clay layer consists of a main aquifer separated from a thin phreatic top layer by an ensemble of fine silty layers. The stratigraphy of the site is given in Table 1.It is noted that the course of the River Aa has been modified in the twentieth century, which explains the absence of alluvial strata corresponding to the River Aa itself.
Based on head data of the Dutch Geological Survey within 5 km of the field site, the groundwater system is a recharge area, drained by the River Aa and its tributary streams. It is a rural area, mainly covered by crop fields and meadows, with occasional patches of woods.
A map of the River Aa and the piezometers is shown in Fig. 1. Heads and stream levels were measured with pressure transducers. Piezometers P7 and P8 were screened at 4 m bgl and are located at a distance of 25 and 50 m from the riverbank, respectively. Piezometers P11 and P12 were screened at 1.5 m bgl and are located at a distance of 25 and 70 m from the riverbank, respectively. The head regularly dropped below the bottom of piezometer P11.
The river stage was recorded 300 m upstream of the piezometers. The precipitation time series was obtained by interpolating the measurements at three weather stations within 15 km from the investigation site. The evaporation time series was obtained from a weather station 11 km from the field site.
The evaporation values correspond to the Makkink reference evaporation, which is representative for Dutch meadowland cover under average meteorological conditions (Bartholomeus et al. 2013). The measurements in the piezometers, the measured rainfall, evaporation and river stage are used to estimate aquifer parameters to be used in a numerical model of the area.

Method
Response function from a one-dimensional (1D) model schematization For application of the flood-wave method in a time series framework, a vertical cross-section is considered along the dashed line in Fig. 1. The cross-section is shown in Fig. 2. The conceptual model consists of a thin phreatic top layer consisting of relatively low permeability material underlain by a semi-pervious layer (aquitard), and a semi-confined layer. The following approximations are adopted.
& The stream fully penetrates the aquifer. Head loss due to stream-line contraction or due to a semi-pervious streambed are lumped in the specific resistance of the streambed (resistance per unit length of streambed) w [TL −1 ] defined as: where Q s [L 2 T −1 ] is the flux from the aquifer to the stream per unit length of stream, h(x = 0) [L] is the head at the interface between the semi-pervious stream bank and the aquifer, and h s [L] is the stream stage. & The boundary opposite to the river is approximated by a zero constant head boundary, at a distance 2 L from the stream. For the case of precipitation and evaporation, this is equivalent to a water divide at a distance L from the stream. The heads in the two layers satisfy the following set of two linked differential equations: where h 1 (t) and h 2 (t) are the heads in the phreatic and semiconfined layers, respectively, x is the distance from the stream bank, R [LT −1 ] is the areal recharge, T[L 2 T −1 ] is the transmissivity of the semi-confined layer, S[−] is the storage coefficient of the phreatic layer, c[T] is the resistance to vertical flow of the aquitard. For the step response to recharge, the boundary conditions are: with h s = 1 for t > 0. For both step responses, the initial conditions are: Solutions for the two step responses are obtained with a Laplace transformation. The Laplace transformation of the differential equation and associated boundary conditions are given in the electronic supplemental material (ESM). The solution in the Laplace domain for the step response to precipitation is and where h 1 and h 2 are the Laplace transformed step responses in layers 1 and 2, p is the Laplace parameter, and γ is The responses to precipitation and evaporation are assumed to be equal in magnitude but opposite in sign.
For the step response to the river stage, and These solutions can be verified by substituting them in the corresponding differential equations and boundary conditions. Back transformation of the step functions from the Laplace domain to the time domain is performed numerically by the method of Stehfest (1970).

Time-series modeling
The extended flood-wave method, now in a time series framework, is run by calculating the groundwater heads at each time step and at each piezometer using Eq. (6). The parameters of the time series model are estimated by a modified Gauss-Newton algorithm (Hill 1998) by maximizing the Nash-Sutcliff coefficient E (Nash and Sutcliffe 1970) defined as: , and parameter α of the exponentially decreasing noise model. These parameters are estimated by maximizing the Nash Sutcliffe coefficient (Eq. 17). The drainage base is fixed to the average stream stage over the simulation period. The river-stage time series was consequently modified by taking the stage relative to the average stage instead of taking the absolute stage value.

Analysis and interpretation
The observed heads and the heads explained by the deterministic part of the time series model are presented in Fig. 3, while the separate contributions of the three stresses are presented in Fig. 4. The observed heads show that the average head and the amplitude of the fluctuation increases with the distance from the draining stream and is the largest for the phreatic piezometer P12 located 70 m from the stream bank (Fig. 1). This is satisfactorily reproduced by the time series model. The peaks observed for P12 caused by precipitation cannot be simulated by the time series model. For the semi-confined piezometers P7 and P8, the influence of the stream stage dominates the fast fluctuations, while precipitation and evaporation cause slower fluctuations. Within the modeled time period (October 2011-February 2012), evaporation decreases which is reflected by the decreasing contribution (in absolute value) of the evaporation component for all three piezometers. Note also that the contribution of the stream stage to the head fluctuations in piezometer P8 is damped compared to piezometer P7 and is almost absent in the case of the phreatic piezometer P12.
The Nash-Sutcliffe coefficients are given for each piezometer in Table 2. The second column gives the coefficients including the noise model, while the third column gives the coefficients for the deterministic part only.The optimal parameters are given in the following. The estimates of the 95 % confidence intervals are given in brackets. The estimated correlation coefficients are given in Table 3.  The confidence intervals vary from ± 21 % to ± 50 %, which is similar to confidence intervals obtained with pumping tests. The distance L is strongly correlated with the noise decay parameter α. Note, however, that the two parameters are not estimated for use in a distributed numerical groundwater model, like the other estimated parameters.

Transmissivity
The value of the transmissivity of the semi-confined layer is lower than a priori expected. The layer thickness as described in Table 1 is estimated as 25 m, and from the bore log descriptions, a hydraulic conductivity of a least 10 m d −1 is expected. An estimated transmissivity of 108 m 2 d −1 suggests an aquifer thickness less than 15 m. An explanation for an apparently thinner aquifer is the presence of silt and peat layers in the formation of Stramproy, constituting the lower 10 m of the semi-confined layer. These semi-pervious layers were assumed to be discontinuous with no confining effect, but the results suggest that the formation of Stramproy acts as a semi-pervious layer reducing the thickness of the investigated semi-confined layer.

Storage coefficient
The value of the storage coefficient obtained for the phreatic layer is 0.14, which is a reasonable value for phreatic layers. It is interesting to mention that Barlow et al. (2000) applied the flood-wave method to find a specific storage coefficient of 9.8 × 10 −5 m −1 for a shallow water-table aquifer with a thickness of about 20 m. The explanation given for this apparent elastic storage was that the thick capillary fringe confines the aquifer. This does not seem to be the case in the present study.
An important difference is that the filter in Barlow et al. (2000) is much deeper than the filter used in this study.

Specific streambed resistance
The response functions lump the head loss due to a semipervious streambed and head loss due to the significant vertical flow component in the vicinity of the stream; the latter is referred to as stream-line contraction. The estimated value of the specific streambed resistance is 0.044 d m −1 . This low value suggests that the head loss is exclusively the result of stream-line contraction. This is supported by field observation of the streambed, which did not reveal the presence of a semipervious riverbed. This hypothesis is tested by evaluating the magnitude of the head loss due to stream-line contraction using an analytical, two-dimensional (2D) cross-sectional model of an aquifer discharging into a stream.
The 2D cross-section shown in Fig. 5 represents an aquifer fed with areal recharge R, discharging into a shallow stream. The aquifer is a finite strip of thickness D [L], with horizontal and vertical hydraulic conductivities k x and k y , respectively. The bottom and top boundaries are impermeable, except along the streambed. The width of the stream is 2B [L]. The origin of the coordinate system is at the stream-aquifer boundary.
The stream is in direct contact with the aquifer, with no semi-pervious streambed. Aquifer discharge into the stream is assumed to be equally distributed over the streambed. The solution for the head φ(x,y) relative to the stream stage is: The derivation of this solution is given in the ESM. The 2D head distribution is plotted in Fig. 6 for an isotropic situation with a recharge rate of R = 0.001 m d −1 , L = 640 m, thickness D = 15 m, and k x = k y = 6.7 m d −1 .
The 2D cross-sectional model is used to estimate the magnitude of the head loss due to stream-line contraction. The 1D (Dupuit) solution to the stated problem with Eq. (9) and h s = 0, is well known and equals Head loss by stream-line contraction is accounted for by the specific streambed resistance w. Different values of w correspond to different values of the vertical anisotropy factor. For example, consider an aquifer of thickness D = 15 m and horizontal permeability k x = 6.7 m d −1 (transmissivity T = 108 m 2 d −1 ). The head calculated by the 1D and 2D models are compared in Fig. 7. Blue corresponds to the head of the 1D model, red corresponds to the head of the isotropic 2D model at a depth of half of the layer thickness. The value of w in the 1D model is adjusted so that the 1D and 2D models coincide for large values of x, which gives w = 0.04 d m −1 . This value is approximately equal to the value obtained by parameter estimation of the time series model so that it may be concluded that the specific streambed resistance estimated with the extended flood-wave method can be used to represent head losses due to stream-line contraction.

Use of the derived parameters in a numerical groundwater model
When a pumping test has been performed and interpreted using an analytical well function, it is a standard practice to enter the transmissivity obtained into a numerical model of the tested aquifer. Similarly, in this section, the aquifer parameters estimated with the extended flood-wave method are used in a distributed numerical groundwater flow model of the investigated field site.
A distributed numerical groundwater model of the field site was built using the same schematization and approximations. The numerical model was implemented with the finite element code MicroFEM (Hemker and de Boer 1997), which allows for the refinement of the mesh along the streams, which was imported from a GIS shape file. The numerical model consists of a phreatic, low-permeability top layer overlying a semiconfined layer where the Dupuit approximation is adopted. Horizontal flow in the phreatic layer is made negligible by fixing the transmissivity to a small value. Model boundaries are either head-dependent when corresponding to a stream, or no-flow boundaries when approximately corresponding to a water divide. Head-dependent boundaries are attributed a head value of zero assuming that stream fluctuations do not influence each other. The modeled area is shown in Fig. 1.
The streambed resistance w′ [T] in the numerical model is obtained through multiplication of the specific resistance w [TL −1 ] obtained with the extended flood-wave method through multiplication with the half width B of the stream with B = 7 m.
Parameters of the numerical model that were not estimated with the extended flood-wave method are: -Streambed resistance for streams other than the River Aa: 0.5 d -Width of streams other than the River Aa: 2 m -Transmissivity of the phreatic layer: 0.1 m 2 d −1 -Specific storage coefficient of the semi-confined layer: 10 −4 m −1 Trying other realistic values for these parameters had only a minor impact on the calculated heads.
Fig. 5 Two-dimensional groundwater flow of precipitation discharging into a stream Groundwater heads are calculated with MicroFEM as suggested by Olsthoorn (2008) by first evaluating the step responses at the place of the piezometers, after which head fluctuations are obtained by convolution of the block response functions with their corresponding stress time series. Finally, the simulated groundwater heads are obtained by adding the drainage base (which is the mean river stage here) as given in Eq. (6). The fit obtained from the numerical groundwater model is satisfying for the semi-confined piezometers, with Nash-Sutcliff coefficients of 0.87 and 0.80 for P7, P8 respectively. The fit for the phreatic piezometer P12 is less good, with a Nash-Sutcliffe coefficient of 0.73, similar to the value of 0.76 obtained with the extended flood-wave method ( Table 2). As with the extended flood-wave method, the model failed to reproduce the fluctuation peaks in piezometer 12 (Fig. 8).

Discussion and conclusion
The objective of this study is to derive aquifer parameters for use in groundwater models from naturally fluctuating heads observed in the vicinity of a stream. The original flood-wave method cannot be applied when the effects of stream stage variations cannot be distinguished from those of precipitation and evaporation by simple inspection of the groundwater head hydrograph. To deal with this problem, the flood-wave method is implemented in the framework of time series analysis to identify the fluctuations associated with each of the stresses (in this paper: precipitation, evaporation, and stream stage variations). The method is called the extended flood-wave method. Convolution of a stress with its corresponding response function provides the effect of that stress on the head. From a time-series modeling perspective, the method proposed is a variation of the method of predefined response functions (Von Asmuth et al. 2002). The response functions of the extended flood-wave method are to be compared with the well function of a pumping test: they translate observed heads into aquifer parameters with a minimum of parameters. An important difference with the original flood-wave method and pumping tests is that aquifer parameters are estimated from the superimposed effects of precipitation, evaporation, and stream stage fluctuations.
The method is illustrated with a case study for an aquifer drained by a lowland river in the Netherlands. The response functions of the time series model represent a cross-section of an aquifer underlying a low-permeability phreatic layer, discharging into a stream. The model describes the essential features of the hydrogeological situation, while keeping it as simple as possible to restrict the number of parameters to optimize. The time series model results in a good fit for the semi-confined piezometers and reproduces the slow fluctuations of the phreatic top layer, but fails to reproduce the quick reactions in the top layer, probably due to nonlinear processes which are not taken into account by the model.
The order of magnitude of the estimated parameters gives qualitative insight into the groundwater system considered. The value of the transmissivity, for example, suggests a new interpretation of the bore logs. The intercalated silt and peat sub layers, revealed by the bore logs at a depth of about 15 m below ground level, might practically form the aquifer bottom instead of a deeper clay layer as initially assumed. The low value found for the resistance of the streambed suggests the absence of a semi-pervious riverbed. Head loss is the result of stream-line contraction in the vicinity of the river, as confirmed by comparing head losses evaluated with an analytical solution for 2D flow in a vertical cross section of an aquifer discharging into a stream.
As for pumping tests, aquifer parameters that are estimated with the extended flood-wave method can be used in a numerical distributed groundwater flow model as prior estimates. It Some evaluative remarks are made about the methodology proposed in this paper. First, the time series model was fitted over a relatively short time period which did not allow the observations time series to be split into a calibration and a validation period. Note that this is similar to pumping tests Second, the conceptual model needs to be kept as simple as possible while incorporating sufficient complexity to match the hydrogeological situation. In an early phase of this study, a simpler groundwater model without the phreatic layer was used, but no reasonable fit with the observed head was possible. The minimum complexity that needs to be incorporated is the additional layer with phreatic storage. Third, the extended flood-wave method relies on a simplification of the reality like any model. The validity of the approximations needs to be considered by the practitioner for each new situation. For example in this study, the fluctuations of the river had negligible impact on the distance between the observation wells and the riverbank. This might not be the case for other rivers. The parallel should be drawn again with pumping tests requiring the choice of an adequate well function. Different contexts require different solutions. Barlow et al. 2000 offer a number of solutions that could be used as an alternative.