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

Abstract

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.

Résumé

La méthode de l’onde de crue est appliquée dans le cadre d’une analyse de séries de niveaux piézométriques afin d’estimer des paramètres de modèles d’écoulement d’eau souterraine. La méthode qui en résulte est applicable dans des situations où les fluctuations de niveaux piézométriques à proximité d’un cours d’eau dépendent aussi de façon significative des régimes de précipitation et d’évaporation. Les fonctions de réponse dans l’analyse de séries de niveaux piézométriques sont générées par un modèle analytique d’écoulement d’eau souterraines décrivant les interactions entre un cours d’eau et un aquifère. Ces fonctions de réponse analytiques jouent le même rôle que les fonctions de réponse dans un essai de pompage, c’est à dire traduire les variations de niveaux piézométriques observées en termes de paramètres de modèle d’écoulement d’eau souterraine au moyen d’une expression analytique parcimonieuse. Une différence importante par rapport à la méthode de l’onde de crue traditionnelle et aux essais de pompage est que les paramètres du modèle d’écoulement d’eau souterraine sont inférés des effets combinés du régime de précipitation, d’évaporation et des fluctuations de niveau d’un cours d’eau. Les fluctuations de niveaux piézométriques naturelles sont séparées en contributions, chacune liée à une cause particulière. La méthode est illustrée en utilisant des niveaux piézométriques mesurés à proximité d’un cours d’eau de basse altitude aux Pays-Bas. Une attention particulière est prêtée à l’interprétation de la résistance à l’écoulement du lit du cours d’eau qui apparait être le résultat des courbures des lignes d’écoulement de l’eau souterraine et non pas de la présence d’un lit de cours d’eau peu perméable. Cette conclusion est tirée en calculant les pertes de hauteur piézométriques au moyen d’un modèle analytique dans un plan vertical en deux dimensions.

Resumen

El método de la onda de crecida se lleva a cabo en el marco del análisis de series de tiempo para estimar parámetros del acuífero para su uso en un modelo de agua subterránea. El método de la onda de crecida ampliada es aplicable a situaciones en que las fluctuaciones del agua subterránea están afectadas significativamente por la precipitación y la evaporación variables con el tiempo. Las funciones de respuesta para el análisis de las series de tiempo se generan con un modelo analítico de agua subterránea que describe la interacción agua superficial-acuífero. Las funciones de respuesta del análisis juegan el mismo papel que la función de pozo en un ensayo de bombeo, que es traducir las variaciones observadas en la carga hidráulica en los parámetros del modelo del agua subterránea por medio de una ecuación parsimoniosa del modelo. Una diferencia importante en comparación entre las pruebas de los métodos tradicionales de onda de crecida y las de bombeo del acuífero es que los parámetros se deducen de los efectos combinados de las fluctuaciones de la precipitación, la evaporación y el nivel de la corriente. Las fluctuaciones que ocurren naturalmente se separan de las contribuciones de diferentes cargas. El método propuesto se ilustra con los datos recogidos cerca de un río de tierras bajas en Holanda. Se pone especial énfasis en la interpretación de la resistencia en el lecho del río. La resistencia del lecho del río es el resultado de la contracción de la línea de la corriente en lugar de la capa semipermeable en el cauce del río, la cual se concluye mediante la comparación con la pérdida de carga hidráulica calculada en una sección transversal de un modelo analítico bidimensional.

摘要

在时间序列分析框架内实施洪波法以估算地下水模型中的含水层参数。在地下水波动受到随时间变化的降水和蒸发显著影响的情况下,扩展的洪波法非常适用。生成了时间序列的响应函数,以及描述河流-含水层相互作用的解析地下水模型。解析响应函数与抽水实验中的井函数发挥着同样的作用,井函数通过简约的模型方程式将观测到的水头变化转换成地下水模型参数。与传统的洪波法和抽水实验相比,一个重要的差别就是含水层参数根据降水量、蒸发量和河流水位波动推断出来。从不同的压力中分离出自然出现的波动所占的比重。根据荷兰一个低地河流附近收集到的资料描述了所提出的方法。特别强调了河床阻力的解译。河床的阻力是河流线收缩造成的,而不是半渗透的河床造成的,这是通过与解析二维横断面模型计算出的水头损失进行比较得出的结论。

Resumo

O método da onda de inundação foi implementado dentro de uma abordagem de series temporais para estimar parâmetros de aquíferos para uso em modelos de águas subterrâneas. O método da onda de inundação estendido resultante é aplicável em situações onde as flutuações das águas subterrâneas são significativamente afetadas pela precipitação e evapotranspiração variante no tempo. Funções de resposta para séries temporais são geradas com um modelo analítico de águas subterrâneas descrevendo a interação rio-aquífero. Funções de resposta analíticas desempenham o mesmo papel que a função do poço em um teste de bombeamento, que é o de traduzir variações de carga piezométrica em parâmetros de modelos de águas subterrâneas pelo uso de uma equação de modelo parcimoniosa. Uma diferença importante quando comparado ao método da onda de inundação tradicional e testes de bombeamento é que os parâmetros do aquífero são inferidos pelo efeito combinado da precipitação, evapotranspiração e flutuações no nível do rio. Flutuações ocorridas naturalmente são separadas em contribuição por diferentes estresses. O método proposto é ilustrado com dados coletados próximos a um polder de rio nos Países Baixos. Ênfase especial é dedicada à interpretação da residência no leito do rio. A residência do leito do rio é resultado de uma contração da linha de fluxo ao invés de um leito de rio semipermeável, conclusão que foi feita através de comparação com a perda de carga calculada com um modelo bidimensional em transepto.

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-series-analysis 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):

$$ \varphi (t)={\displaystyle \underset{0}{\overset{t}{\int }}p\left(\tau \right)\theta \left(t-\tau \right)d\tau } $$

(1)

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

$$ s(t)={\displaystyle \underset{0}{\overset{t}{\int }}\theta \left(t-\tau \right)d\tau } $$

(2)

The step function has the dimension of length per dimension of stress.

The function

$$ \psi \left(t,\Delta t\right)=s(t)-s\left(t-\Delta t\right) $$

(3)

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:

$$ \varphi \left({t}_j\right)={\displaystyle \sum_{i=1}^j{p}_i}\psi \left({t}_j-{t}_{i-1},\Delta {t}_i\right) $$

(4)

where

$$ {t}_j={\displaystyle \sum_{i=1}^j\Delta {t}_i} $$

(5)

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\left({t}_j\right)=d+{\varphi}_{\mathrm{p}}(t)+{\varphi}_{\mathrm{e}}(t)+{\varphi}_{\mathrm{s}}(t)+n(t) $$

(6)

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, 67-km-long lowland river, with an average flow of 11 m³s⁻¹ at its mouth.

Fig. 1

Location of the field site at different scales: in the Netherlands, in respect to local streams, in respect to the River Aa

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.

Table 1 Site stratigraphy

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.

Fig. 2

Conceptual model of the vertical cross-section for response to recharge

• 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⁻¹] defined as:

  $$ {Q}_{\mathrm{s}}=\frac{h\left(x=0\right)-{h}_{\mathrm{s}}}{w} $$

  (7)

  where Q _s [L²T⁻¹] 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 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).

• 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.

• Flow in the top layer is vertical.

• The resistance to vertical flow is neglected in the semi-confined 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:

$$ \left\{\begin{array}{l}\frac{h_2-{h}_1}{c}=S\frac{\partial {h}_1}{\partial t}-R\\ {}T\frac{\partial^2{h}_2}{\partial {x}^2}+\frac{h_1-{h}_2}{c}=0\end{array}\right. $$

(8)

where h ₁(t) and h ₂(t) are the heads in the phreatic and semi-confined layers, respectively, x is the distance from the stream bank, R [LT⁻¹] is the areal recharge, T[L²T⁻¹] 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:

$$ \begin{array}{l}\frac{\partial {h}_2}{\partial x}=0,\ x=L\\ {}T\frac{\partial {h}_2}{\partial x}=\frac{h_2-{h}_{\mathrm{s}}}{w},\ x=0\end{array} $$

(9)

where the stream stage h _s [L] is 0 m.

For the step response to stream stage fluctuations, the boundary conditions are:

$$ \begin{array}{l}{h}_2=0,\ x=2L\\ {}T\frac{\partial {h}_2}{\partial x}=\frac{h_2-{h}_{\mathrm{s}}}{w},\ x=0\end{array} $$

(10)

with h _s = 1 for t > 0.

For both step responses, the initial conditions are:

$$ {h}_1\left(x,t=0\right)={h}_2\left(x,t=0\right)=0 $$

(11)

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

$$ {\overline{h}}_1=\frac{p{\overline{h}}_2+c}{p\left(cSp+1\right)} $$

(12)

and

$$ {\overline{h}}_2=\frac{1}{Sp}\left(\frac{- \cosh \left(\gamma \left(L-x\right)\right)}{p\left(Tw\gamma \sinh \left(\gamma L\right)+ \cosh \left(\gamma L\right)\right)}+\frac{1}{p}\right) $$

(13)

where ${\overline{h}}_1$ and ${\overline{h}}_2$ are the Laplace transformed step responses in layers 1 and 2, p is the Laplace parameter, and γ is

$$ \gamma =\sqrt{\frac{Sp}{T\left(cSp+1\right)}} $$

(14)

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,

$$ {\overline{h}}_1=\frac{{\overline{h}}_2}{cSp+1} $$

(15)

and

$$ {\overline{h}}_2=-\frac{ \sinh \left(\gamma \left(2L-x\right)\right)}{p\left(Tw\gamma \cosh \left(2\gamma L\right)+\left( \sinh \left(2\gamma L\right)\right)\right)} $$

(16)

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:

$$ \operatorname{E}=1-\frac{{\displaystyle \sum_{i=1}^N{\left({h}_{\mathrm{o},i}-{h}_{\mathrm{m},i}\right)}^2}}{{\displaystyle \sum_{i=1}^N{\left({h}_{\mathrm{o},i}-{\mu}_{\mathrm{o}}\right)}^2}} $$

(17)

where h _o,i is the observed head at time i, h _m,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 [L²T⁻¹], 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⁻¹], 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.

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.

Fig. 3

Time series model of the head at piezometers P7, P8, and P12 with observed head, observed river stage, and observed precipitation

Fig. 4

Modeled contributions of the different stresses to the fluctuations of the groundwater head at piezometers P7, P8, and P12, as inferred from the time series model

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.

Table 2 Nash-Sutcliffe coefficient of the modeled fluctuations using the extended flood-wave method

Table 3 Correlation coefficients for the parameters of the extended flood-wave method

• Transmissivity semi-confined layer T : 108 m² d⁻¹ [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⁻¹ [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.

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⁻¹ is expected. An estimated transmissivity of 108 m² d⁻¹ 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⁻⁵ m⁻¹ 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 semi-pervious 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⁻¹. 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 semi-pervious 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.

Fig. 5

Two-dimensional groundwater flow of precipitation discharging into a stream

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:

$$ \varphi \left(x,y\right)={\displaystyle \sum_{n=1}^{\infty}\frac{\frac{d{\overline{\varphi}}_n}{dy}\left(y=0\right)}{\alpha n \sinh \left(\alpha nD\right)}\left( \cosh \left(\alpha n\left(D+y\right)\right) \cos \left(\frac{\alpha n\left(L-B-x\right)}{\sqrt{a}}\right)- \cos \left(\frac{\alpha n}{\sqrt{a}}\left(1-B/L\right)\right) \cosh \left(\alpha nD\right)\right)} $$

(18)

with $a=\frac{k_x}{k_y}$, $\alpha n=\frac{n\pi }{L}\sqrt{a}$ and $\frac{\mathrm{d}{\overline{\varphi}}_n}{\mathrm{d}y}\left(y=0\right)=-2{\left(-1\right)}^n\frac{R}{ky}\frac{\sqrt{a}}{\alpha nB} \sin \left(\frac{\alpha nB}{\sqrt{a}}\right)$.

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⁻¹, L = 640 m, thickness D = 15 m, and k _x = k _y = 6.7 m d⁻¹.

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⁻¹, L = 640 m, thickness D = 15 m, and k _x  = k _y  = 6.7 m d⁻¹; the contour line h = 0 is the approximate bottom of the stream

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

$$ \varphi (x)=R\left(-\frac{x^2}{2{k}_xD}+\frac{Lx}{k_xD}+Lw\right) $$

(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 k _x = 6.7 m d⁻¹ (transmissivity T = 108 m² d⁻¹). 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⁻¹. 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.

Fig. 7

Heads for 1D model with w = 0.04 d m⁻¹ and recharge R = 0.001 m d⁻¹ (blue) and 2D model (red) at depth of half of layer thickness

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 semi-confined 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⁻¹] 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² d⁻¹

• Specific storage coefficient of the semi-confined layer: 10⁻⁴ m⁻¹

Trying other realistic values for these parameters had only a minor impact on the calculated heads.

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. 8

Comparison of heads simulated with the time series model and heads calculated with the numerical groundwater flow model with the same parameters as the time series model

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.