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

Hydrogeology Journal, Jun 2016

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs10040-016-1436-5.pdf

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

A time-series analysis framework for the flood-wave method to estimate groundwater model parameters Christophe Obergfell 0 1 2 Mark Bakker 0 1 2 Kees Maas 0 1 2 0 KWR Watercycle Research Institute , Nieuwegein , The Netherlands 1 Water Resources Section, Delft University of Technology , Delft , The Netherlands 2 Christophe Obergfell 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 timeseries 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. The Netherlands; Time series analysis; Groundwater/surface-water relations; Analytical solutions; Numerical modeling 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) . 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): φð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 t Z 0 sðtÞ ¼ θðt−τ Þdτ The step function has the dimension of length per dimension of stress. The function ψðt; ΔtÞ ¼ sðtÞ−sðt−ΔtÞ 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 Δti is the length of stress period i. Stress pi is constant over stress period i from t = ti to t = ti + Δti. Since the system is approximated as linear, the head at time tj can be obtained by summing the effects at time tj of all past stress periods: j φ t j ¼ X piψ t j−ti−1; Δti i¼1 where j t j ¼ X Δti i¼1 For application of the extended flood-wave method in this paper, the heads h(t) in the aquifer are obtained as the following sum: h t j ¼ d þ φpðtÞ þ φeðtÞ þ φsðtÞ þ nðtÞ 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 ð2Þ ð3Þ ð4Þ ð5Þ ð6Þ 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 m3s−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 North Sea Germany Netherlands 0 40 80 km Belgium Catchment area river Aa Legend Isohypse of estimated mean head in semiconfined layer (m above sea level) River Aa Secondary stream Field site Cross-section of 1D groundwater flow model Piezometer Pdroirjeeccttiioonn B A Distance L 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: Qs ¼ hðx ¼ 0Þ−hs w ð7Þ & & & where Qs [L2T−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 hs [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 piezometers are approximately positioned along a flow line. Precipitation surplus reaches the water table instantaneously (the depth to the water table is about 1 m). a bgl below ground level, with ground level at 2.5 m NAP (NAP is the Dutch datum, approximately corresponding to mean sea level) b Base of groundwater model & & & The base of the system is impermeable. The storage of the semi-confined layer is negligible with respect to the phreatic storage of the phreatic layer. The semi-confined layer has a uniform transmissivity. Fig. 2 Conceptual model of the vertical cross-section for response to recharge z Stream (River Aa) hs Legend hs x Semi-confined layer Impermeable base Stream stage & & & Flow in the top layer is vertical. The resistance to vertical flow is neglected in the semiconfined layer (Dupuit approximation). The river stage variations result in negligible changes in the distance between the riverbank and the piezometers. The heads in the two layers satisfy the following set of two linked differential equations: 8 h2−h1 <> c ∂2h2 >: T ∂x2 þ ∂h1 ¼ S ∂t h1−h2 c −R where h1(t) and h2(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[L2T−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: ∂h2 ∂x ¼ 0; x ¼ L T ∂∂hx2 ¼ h2w−hs ; x ¼ 0 where the stream stage hs [L] is 0 m. h2(x,t) h1(x,t) h1(x,t) Low -permeability phreatic layer Aquitard Head in phreatic layer h2(x,t) Semi-pervious Stream bank Initial stream stage and initial heads Head in semiconfined layer Water divide L ð8Þ ð9Þ For the step response to stream stage fluctuations, the boundary conditions are: h2 ¼ 0; x ¼ 2L T ∂∂hx2 ¼ h2w−hs ; x ¼ 0 with hs = 1 for t > 0. For both step responses, the initial conditions are: h1ðx; t ¼ 0Þ ¼ h2ðx; t ¼ 0Þ ¼ 0 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 ð10Þ ð11Þ ð12Þ ð13Þ ð15Þ ð16Þ ph2 þ c h1 ¼ pðcSp þ 1Þ and γ ¼ sffiffiffiffiffiffiffiffiffiSffiffiffipffiffiffiffiffiffiffiffiffiffiffi T ðcSp þ 1Þ 1 −coshðγðL−xÞÞ 1 h2 ¼ Sp pðT wγsinhðγLÞ þ coshðγLÞÞ þ p where h1 and h2 are the Laplace transformed step responses in layers 1 and 2, p is the Laplace parameter, and γ is ð14Þ Analysis and interpretation 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, h2 h1 ¼ cSp þ 1 and h2 ¼ − sinhðγð2L−xÞÞ pðT wγcoshð2γLÞ þ ðsinhð2γLÞÞÞ 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 GaussNewton algorithm (Hill 1998) by maximizing the NashSutcliff coefficient E (Nash and Sutcliffe 1970) defined as: E ¼ 1− i¼1 N X N X i¼1 ho;i−hm;i ho;i−μo 2 2 ð17Þ where ho,i is the observed head at time i, hm,i is the modeled head at time i, and μo is the average observed head. It is recalled that the parameters of the extended flood-wave method are the transmissivity of the semi-confined layer T [L2T−1], the storage coefficient of the phreatic layer S [−], the resistance to vertical flow of the aquitard c [T], the specific resistance of the streambed w [TL−1], the distance between the riverbank and the constant head boundary 2L [L] and the drainage base d [L], 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. 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. – Transmissivity semi-confined layer T : 108 m2 d−1 [80–147] – – – – – Resistance to vertical flow of the aquitard c : 79 d [48–127] Phreatic storage coefficient S : 0.14 [0.11–0.17] Streambed specific resistance w : 0.044 d m−1 [0.031–0.065] Distance L : 640 m [420–986] Exponent of noise model α : 0.15 [0.11,0.19] 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. Table 2 Nash-Sutcliffe coefficient of the modeled fluctuations using the extended flood-wave method 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 m2 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. ∞ φðx; yÞ ¼ X n¼1 αnsinhðαnDÞ ddφyn ðy ¼ 0Þ An important difference is that the filter in Barlow et al. (2000) is much deeper than the filter used in this study. Specif ic 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 kx and ky, 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: coshðαnðD þ yÞÞcos αnðL−B−xÞ pffiaffi −cos αn pffiaffi 1−B.L coshðαnDÞ ð18Þ w i t h a ¼ kkyx, nπ pffiaffiffi a n d αn ¼ L ddφyn ðy ¼ 0Þ ¼ −2ð−1Þn kRy αpnffiaBffisin αpnaB . ffiffi 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 kx = ky = 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 hs = 0, is well known and equals x2 Lx φðxÞ ¼ R − 2kxD þ kxD þ Lw ð19Þ 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 kx = 6.7 m d−1 (transmissivity T = 108 m2 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 Fig. 5 Two-dimensional groundwater flow of precipitation discharging into a stream 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 m2 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. 0.0 -5.0 ) m ( h t p e D -10.0 -15.0 -5.0 0.0 5.0 distance to stream bank (m) Fig. 6 Groundwater head (in meters) in a 2D cross-section of an aquifer discharging into a stream, obtained from Eq. (18) with a recharge rate of R = 0.001 m d−1, L = 640 m, thickness D = 15 m, and kx = ky = 6.7 m d−1; the contour line h = 0 is the approximate bottom of the 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). Fig. 7 Heads for 1D model with w = 0.04 d m−1 and recharge R = 0.001 m d−1 (blue) and 2D model (red) at depth of half of layer thickness 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 is essential that the numerical model shares the same schematization and assumptions as used in the extended flood-wave method, similar to what is done with pumping tests. A numerical groundwater model, parameterized in this way, results in a good fit, except again for the quick reactions in the top layer. 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 that are usually conducted over a short period of time. A validation period is particularly recommended when a time series model is used for predictions. 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. Acknowledgements The authors thank Artesia Consultants, (Schoonhoven, The Netherlands) and the Water Board Aa & Maas (‘sHertogenbosch, The Netherlands) for providing all field data and for their constructive comments. This work was performed in the cooperation framework of Wetsus, European centre of excellence for sustainable water technology (www.wetsus.nl). Wetsus is co-funded by the Dutch Ministry of Economic Affairs and Ministry of Infrastructure and Environment, the European Union Regional Development Fund, the Province of Fryslân and the Northern Netherlands Provinces. The authors would like to thank the participants of the research theme ‘Groundwater technology’ for the fruitful discussions and their financial support. Open Access This article is distributed under the terms of the Creative C o m m o n s A t t r i b u t i o n 4 . 0 I n t e r n a t i o n a l L i c e n s e ( h t t p : / / creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Barlow PM , DeSimone LA , Moench AF ( 2000 ) Aquifer response to stream-stage and recharge variations, II: convolution method and applications . J Hydrol 230 : 211 - 229 . doi: 10 .1016/S0022- 1694 ( 00 ) 00176 - 1 Bartholomeus RP , Voortman BR , Witte JPM ( 2013 ) Measurements and knowledge of hydrological processes are required for accurate estimations of evaporation with groundwater models (in Dutch) . Stromingen 19 ( 2 ): 37 - 52 Besbes M , de Marsily G ( 1984 ) From infiltration to recharge: use of a parametric transfer function . J Hydrol 74 : 271 - 293 Box GEP , Jenkins GM ( 1976 ) Time series analysis: forecasting and control . Prentice Hall , Englewood Cliffs, NJ Cooper HH , Rorabaugh MI ( 1963 ) Ground-water movements and bank storage due to flood stages in surface streams . US Geol Surv Water Suppl Pap 1536-J Ha K , Koh D-C , Yum B-W, Lee K-K ( 2007 ) Estimation of layered aquifer diffusivity and river resistance using flood wave response model . J Hydrol 337 : 284 - 293 . doi: 10 .1016/j.jhydrol. 2007 . 01 .040 Hall FR , Moench AF ( 1972 ) Application of the convolution equation to stream-aquifer relationships . Water Resour Res 8 : 487 - 493 . doi: 10 .1029/WR008i002p00487 Hemker CJ , de Boer RG ( 1997 ) MicroFEM, 4 . 10 . 42 finite element method groundwater flow simulations computer program . MicroFEM , Amsterdam Hespanha JP ( 2009 ) Linear systems theory . Princeton University Press, Princeton, NJ Hill MC ( 1998 ) Methods and guidelines for effective model calibration . US Geol Surv Water Resour Invest Rep 98-400 Knight JH , Rassam DW ( 2007 ) Groundwater head responses due to random stream stage fluctuations using basis splines . Water Resour Res 43. doi:10.1029/2006WR005155 Moench AF , Barlow PM ( 2000 ) Aquifer response to stream-stage and recharge variations, I: analytical step-response functions . J Hydrol 230 : 192 - 210 . doi: 10 .1016/S0022- 1694 ( 00 ) 00175 -X Nash JE , Sutcliffe JV ( 1970 ) River flow forecasting through conceptual models, part I: a discussion of principles . J Hydrol 10 : 282 - 290 . doi: 10 .1016/ 0022 - 1694 ( 70 ) 90255 - 6 Obergfell C , Bakker M , Zaadnoordijk W , Maas K ( 2013 ) Deriving hydrogeological parameters through time series analysis of groundwater head fluctuations around well fields . Hydrogeol J 21 : 987 - 999 . doi: 10 .1007/s10040-013-0973-4 Olsthoorn TN ( 2008 ) Do a bit more with convolution . Ground Water 46 : 13 - 22 Peterson TJ , Western AW ( 2014 ) Nonlinear time-series modeling of unconfined groundwater head . Water Resour Res 50 : 8330 - 8355 . doi: 10 .1002/2013WR014800 Pinder GF , Bredehoeft JD , Cooper HH Jr ( 1969 ) Determination of aquifer diffusivity from aquifer response to fluctuations in river stage . Water Resour Res 5 : 850 - 855 . doi: 10 .1029/WR005i004p00850 Rassam DW , Pagendam D , Hunter H ( 2008 ) Conceptualisation and application of models for groundwater-surface water interactions and nitrate attenuation potential in riparian zones . Environ Model Softw 23 : 859 - 875 Shapoori V , Peterson TJ , Western AW , Costelloe JF ( 2015a ) Decomposing groundwater head variations into meteorological and pumping components: a synthetic study . Hydrogeol J 23 : 1431 - 1448 . doi: 10 .1007 /s10040-015-1269-7 Shapoori V , Peterson TJ , Western AW , Costelloe JF ( 2015b ) Estimating aquifer properties using groundwater hydrograph modelling . Hydrol Process. doi:10 .1002/hyp.10583 Shapoori V , Peterson TJ , Western AW , Costelloe JF ( 2015c ) Top-down groundwater hydrograph time-series modeling for climate-pumping decomposition . Hydrogeol J 23 : 819 - 836 . doi: 10 .1007/s10040-014- 1223-0 Stehfest H (1970) Algorithm 368: numerical inversion of Laplace transforms [D5] . Commun ACM 13 : 47 - 49 . doi: 10 .1145/361953.361969 Venetis C ( 1970 ) Finite aquifers: characteristic responses and applications . J Hydrol 12 : 53 - 62 Von Asmuth JR , Bierkens MFP ( 2005 ) Modeling irregularly spaced residual series as a continuous stochastic process . Water Resour Res 41 , W12404 . doi: 10 .1029/2004wr003726 Von Asmuth JR , Maas K , Bierkens MFP ( 2002 ) Transfer function-noise modeling in continuous time using predefined impulse response functions . Water Resour Res 38 : 1287 . doi: 10 .1029/2001wr001136 Von Asmuth JR , Maas K , Bakker M , Petersen J ( 2008 ) Modeling time series of ground water head fluctuations subjected to multiple stresses . Ground Water 46 : 30 - 40


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs10040-016-1436-5.pdf

Christophe Obergfell, Mark Bakker, Kees Maas. A time-series analysis framework for the flood-wave method to estimate groundwater model parameters, Hydrogeology Journal, 2016, 1807-1819, DOI: 10.1007/s10040-016-1436-5