Parameter estimation for leaky aquifers using the extended
Kalman filter, and considering model and data
measurement uncertainties
H.D. Yeh*, Y.C. Huang
Institute of Environmental Engineering, National Chiao-Tung University, 75 Po-Ai Street, Hsinchu 30039, Taiwan
Received 5 November 2002; revised 22 June 2004; accepted 25 June 2004
Abstract
A method using the extended Kalman filter (EKF) is proposed to identify the hydraulic parameters in leaky aquifer systems both with and without considering the aquitard storage. In the case without considering the aquitard storage, Hantush and Jacob’s model combined with EKF can optimally determine the parameters for the leaky aquifer when analyzing the drawdown data. Coupled with Neuman and Witherspoon’s model, the EKF is also employed to estimate the four parameters of aquifers. The observed drawdown data may be either interpolated using the Lagrangian polynomial or recursively used while implementing the EKF. The proposed method can identify the parameters, using part of the interpolated drawdown data or recursively used data, and obtains results with good accuracy. In the field-pumping test, a long pumping time may not be necessary if the proposed method is implemented on a computer, which is connected to pressure transducers and a data logger. In the process of parameter estimation, the leakage coefficient changes marginally for the first few observations. This phenomenon reflects the fact that there is a time lag between the start of pumping and the leakage effect on the drawdown. The analyses of the data uncertainty demonstrate that the EKF approach is applicable for drawdown data even when it contains white noise or temporal correlated noise. Finally, the choice between Hantush and Jacob’s model and Neuman and Witherspoon’s model depends on the hydrogeological condition of the aquifer system indicated in the analyses of the model uncertainty. Hantush and Jacob’s model is shown to be a good choice for representing the leaky aquifer system if the aquitard storage is comparatively small.
q2004 Elsevier B.V. All rights reserved.
Keywords: Parameter estimation; Kalman filter; Lagrangian polynomial; Groundwater; Leaky aquifer; Model uncertainty
1. Introduction
Hydrogeologic parameters are very important in site characterization, so groundwater hydrologists
often conduct pumping tests to determine hydrogeo-logic parameters, such as hydraulic conductivity and storage coefficient. These parameters are necessary information for quantitative and/or qualitative groundwater studies. In a leaky aquifer the semi-pervious bed (also known as the aquitard), although of very low permeability, may yield significant amounts of water to the adjacent pumped aquifer. As time
www.elsevier.com/locate/jhydrol
0022-1694/$ - see front matter q 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.jhydrol.2004.06.035
* Corresponding author. Fax: C886-35-726050. E-mail address: [email protected] (H.D. Yeh).
increased, leakage across the semi-impervious bed may become appreciable and flow is not restricted to the pumped aquifer alone. The additional water may be derived from storage of the aquitard and adjacent unpumped aquifers. Therefore, the leaky aquifer must be viewed as part of a complex multiple aquifer system.
Two approaches have been developed for dealing with leaky aquifers, one based on the assumption that the aquitard storage is negligible, and the other considering the aquitard storage. Jacob (1946)
developed a partial differential equation for non-steady radial flow in a leaky aquifer, assuming that hydraulic head in the unpumped adjacent aquifer is constant and the storage capacity of the aquitard is negligible. These two assumptions greatly facilitated mathematical treatment of the problem, and Jacob used this approach to develop a large number of solutions to various problems involving flow in aquifers with vertical leakage.
Hantush and Jacob (1955) described non-steady radial flow to a well in a fully penetrated leaky aquifer under a constant pumping rate. In their model, the aquitard is overlain by an unconfined aquifer, and the main aquifer is underlain by an impermeable bed. Their solution is herein called the three-parameter model. Hantush (1960) also presented a modified approach to include the effect of the aquitard storage.
Neuman and Witherspoon (1969) gave a solution describing the drawdown of the lower and pumped aquifer in a hydrogeologic system, which is composed of two confined aquifers and one aquitard. Their solution, which considers the effect of aquitard storage and neglects the drawdown in the unpumped aquifer, is called the four-parameter model. Both the three-parameter and four-parameter models are also mentioned in several recent textbooks, for example,
Dawson and Istok (1991) and Batu (1998).
In the three-parameter model, the graphical method based on Hantush’s or Walton’s type curves (Batu, 1998) requires data plotting work and individ-ual judgment during the curve fitting procedure. Therefore, errors may be introduced during the fitting process. In the four-parameter model, the use of the graphical matching method based on the Neuman and Witherspoon’s model is practically impossible since there will be several families of type curves.
Yeh (1987) used the non-linear least-squares
and finite-difference Newton’s method (NLN) for identifying the parameters of the confined aquifer, a method which has the advantages of high accuracy and quick convergence. Yeh and Han (1989) sub-sequently used NLN to determine the hydraulic parameters of the leaky aquifers.
The Kalman filter was developed by R.E. Kalman in the late 1950s, and its main applications have been in control systems, tracking and navigation of all sorts of vehicles, as well as predictive design of estimation and control systems. Works using the Kalman filter for hydraulic-parameter and water table-related estimations may be divided into two categories. One applies the Kalman filter in a linear system (e.g.Van Geer and Van Der Kloet, 1985; Van Geer and te Stroet, 1990; Van Geer et al., 1990; Lee et al., 2000; Bierkens et al., 2001) and the other deals with non-linear problems using the extended Kalman filter (EKF) (e.g. Chander et al., 1981; Katul et al., 1993; Bierkens, 1998; Cahill et al., 1999).
Chander et al. (1981)estimated the parameters for both non-leaky and leaky aquifers by the iterated EKF. For the leaky aquifer, the measurement equation was a truncated form of Hantush and Jacob’s well function (1955), which was suitable for small leakage coefficient and/or large pumping time. Van Geer and Van Der Kloet (1985) presented two linear, filter-based schemes for parameter estimation in groundwater flow problems. An optimal estimate was simultaneously computed for the original state, i.e., heads and the parameter state.Van Geer and te Stroet (1990)combined MODFLOW with a filtering framework and updated the prior estimates of hydraulic parameter values using an off-line pro-cedure, when minimizing the difference between the actual head measurements and those predicted from the MODFLOW–Kalman filter framework.Van Geer et al. (1990)also adopted the idea of using a filter for state estimation in the absence of significant dynamic behavior and studied the applicability of the filter to a relatively quickly reacting groundwater system.
Katul et al. (1993) applied the EKF to test the determination of the hydraulic conductivity function from a field drainage experiment. Bierkens (1998)
embedded the stochastic differential equation in the EKF algorithm to calibrate the parameters and noise statistics of the stochastic differential equation on a time series of water table depths.Eigbe et al. (1998)
reviewed the problems of using a filter with groundwater flow models, and identified procedures that would improve its efficiency and convenience for field applications.Cahill et al. (1999)derived optimal parameters for an effective large-scale hydraulic conductivity field that considers both the temporal and spatial variations of the moisture content. Lee et al. (2000)presented a linear model to automatically decide the terminal time of a pumping test and to optimally identify the storage coefficient and aniso-tropic aquifer transmissivities and of a confined aquifer, using Cooper and Jacob’s equation and Papadopulos’ approach. Bierkens et al. (2001)
modeled the spatiotemporal variation of shallow water table depth with a regionalized version of an autoregressive exogenous (ARX) time series model. There, the regionalized ARX parameters were estimated by embedding the regionalized ARX model in a space–time Kalman filter.
Beck (1987)reviewed the role of uncertainty in the identification of mathematical models of water quality and in the application of these models to problem of prediction.Eisenberg et al. (1989)selected five broad classes of sources of uncertainty in the area of the geologic repository for radioactive waste. The five types of uncertainty they considered are: (1) systema-tic and random error in measurement; (2) spatial variations in geologic parameters; (3) conceptual model uncertainties relative to geometrical configur-ation, major features, and boundary conditions; (4) physicochemical process modeling; and (5) future states of nature. Freeze et al. (1992) presented an overview and examples of the concepts and tech-niques for dealing with uncertainty in the framework of hydrogeological decision analysis. Wheater et al. (2000) reviewed the problems of uncertainty in the definition of aquifer properties and the scale-dependence of dispersion with specific regard to well-capture zones. Christiaens and Feyen (2001)
compared and evaluated the uncertainties resulting from four ways to obtain soil hydraulic parameters with respect to their resulting uncertainties on different hydrological model outputs. Zheng and Bennett (2002) listed three methods: sensitivity analysis, Monte Carlo method, and the first-order error analysis for evaluating uncertainty.
In this current study, the EKF is used to determine the hydraulic parameters of the leaky aquifer system.
The state vector of the EKF is comprised of the hydraulic parameters and a model, i.e. either Hantush and Jacob’s model or Neuman and Witherspoon’s model, is used as the measurement equation of the EKF. The unknown parameters can be estimated on-line as the measurement data come in and are interpolated by the Lagrangian polynomial. Thus, this proposed approach can be implemented on a computer, which can be connected with a data logger and a pressure transducer that measures the water level in the observation well during the pumping test. The desired hydraulic parameters can thus be estimated in the field on-line and once stable estimates of the parameters are obtained, then the pumping test may be terminated.
The five objectives of this paper are: (1) to develop a numerical approach based on the EKF method which can identify the hydraulic parameters of leaky aquifer in the field on-line; (2) to examine the accuracy and applicability ofChander et al.’s approach (1981)while using the simplified formula of the Hantush and Jacob’s model (1955); (3) to determine four hydraulic parameters using the EKF method based onNeuman and Witherspoon’s model (1969); (4) to test the applicability of the EKF approach if the drawdown data contain white noise or temporal correlated noise; (5) to compare the choice between either Hantush and Jacob’s model or Neuman and Witherspoon’s model for parameter estimation. The data measurement error (herein treated as noise) is related to the problem of measurement data uncertainty, and the choice among various mathematical models for representing a target aquifer belongs to the issue of model uncertainty.
2. Methodology
This section includes two parts: one is the basic framework of the EKF and the other is the Lagrangian polynomial. To implement EKF, the latter approach is applied to interpolate the drawdown data with non-uniform time intervals having a small and uniform time interval.
2.1. Discrete extended Kalman filter
Consider a non-linear dynamical system whose state vector is described as follows
(Grewal and Andrews, 1993)
xkZ f ðxkK1; k K1Þ Cwk (2.1) where xk is the state vector of system at time step
k, f(xkK1,kK1) is the non-linear function of the
system, and wk is the state noise assumed to be
normally distributed with zero mean white (uncor-related) sequence with known covariance Qk. Note
that Qk is assumed to be a zero vector throughout
the filtering steps.
For computing the predicted state estimate, the non-linear implementation equations for the state vector can be described as
^xkðKÞ Z f ð^xkK1ðCÞ; k K1Þ (2.2) where ^xkðKÞ denotes the a priori estimate at k step and ^xkðCÞ epresents the a posteriori estimate at kK1 step. A measurement model of system can be described as (Grewal and Andrews, 1993)
zkZ hðxk; kÞ Cyk (2.3)
where h(xk,k) is the function for the measurement
system and zkis the measurement vector at time step k.
The measurement noise ykis assumed to be a white
sequence with known covariance structure Rk. Note
that the covariance matrix Rk is assumed constant
throughout the filtering process.
As in Eq. (2.2) for computing the predicted measurement, the non-linear implementation equation for the measurement is
^zkZ hð ^xkðKÞ; kÞ (2.4)
where ^zkis predicted measurement vector.
This method requires an initial estimate of the process at some point in time step k, and that estimate is based on the knowledge about the process. The error covariance matrix associated with ^xkðKÞ is also assumed to be known. The prior estimation error ek(K)is defined as xkK^xkðKÞ: Initially, the a priori error covariance matrix of xkK^xkðKÞ; denoted as Pk(K), is
PkðKÞ Z E½ekðKÞe T kðKÞ
Z E b ðxkK^xkðKÞÞðxkK^xkðKÞÞT c (2.5) With the assumption of the a priori error estimate ^xkðKÞ; the measurement zk is used to improve
the a priori estimate, as follows
^xkðCÞ Z ^xkðKÞ C KkðzkK^zkÞ (2.6) where Kkis defined as the Kalman gain and ^xkðCÞ is the updated estimate at step k.
The Kalman filter uses minimum mean-square error as the performance criterion to find the particular
Kk that yields an updated estimate. Therefore, the particular Kalman gain that minimizes the mean-square estimation error is
KkZ PkðKÞH T k½HkPkðKÞH T k CRk K1 (2.7) where the measurement matrix Hkis approximated by
the derivative of estimated measurement ^zk from Eq. (2.4) with respect to each of the state vector. It can be expressed as Hkz vhðx; kÞ vx xZ^xkðKÞ (2.8)
Therefore, Kk can be estimated with known Pk(K),
Hk, and Rk. Detailed derivations of the minimization
process may be found inGrewal and Andrews (1993). The error covariance matrix derived from Eq. (2.5) can also be expressed as
PkðKÞ Z FkK1PkK1ðCÞF T
kK1CQkK1 (2.9)
where FkK1is the state transition matrix. This may be
expressed by the derivative of the state vector estimation equation from Eq. (2.2) as
FkK1z
vf ðx; k K1Þ
vx xZ^xkK1ðKÞ
(2.10)
The error covariance matrix associated with the updated (a posteriori) estimate can be written as PkðCÞ Z fI K KkHkgPkðKÞ (2.11) where Pk(C) denotes the error covariance updated by
the Kalman gain and the prior error covariance and I represents an identity matrix.
Using Eqs. (2.2), (2.4), (2.6), (2.7), (2.9), and (2.11), the recursive process of EKF is then established.
2.2. Lagrangian polynomial
Data where the x-values are not evenly spaced often occur as the result of experimental observations
or when historical data are examined. The Lagrangian polynomial is perhaps the simplest way used to interpolate the unevenly spaced data.
The general form of Lagrangian polynomial is (Gerald and Wheatley, 1994)
PnðxÞ Z XnC1 iZ1 XnC1 i Z 1 i sj x K xi xiKxj f ðxiÞ (2.12)
It is easy to see that the Lagrangian polynomial passes through each of the points used in its construction. Note that this study utilizes a second-degree Lagrangian polynomial to generate interp-olated data points.
3. Application of discrete EKF
This section illustrates how the Kalman filter is coupled with three- and four-parameter models to identify the hydraulic parameters. The difference between the present study andChander et al. (1981)
on the three-parameter model is clearly demonstrated. The three parameter values, i.e. transmissivity T, storage coefficient S, and leakage coefficient L are optimally determined when using the three-parameter model coupled with EKF to analyze the measurement drawdown data. Likewise, the four-parameter model, which has an additional parameter j accounting for the effect of the aquitard storage, combined with EKF can estimate four hydraulic parameters in a similar procedure.
3.1. Leaky aquifer without storage effect in an aquitard
3.1.1. Hantush and Jacob’s model
With reasonable initial guess values for T, S, L, PkK1(C), and Rk, the algorithm of EKF combined
with the Hantush and Jacob’s model can determine the best-fit T, S, and L at each identification step. In their model, the aquitard, with a negligible storage, is overlain by an unconfined aquifer, and the main aquifer is underlain by an impermeable bed. The model of Hantush and Jacob describing the drawdown cone within a leaky aquifer in response to
the pumping as a function of radial distance and time is (Hantush and Jacob, 1955, p.98)
s Z Q 4pTW u; r B (3.1) where s is drawdown, r is the distance between pumping well and observation well, u is dimension-less variable defined as r2S/4Tt, K0 is the vertical conductivity of the aquitard, b0 is thickness of the aquitard, r/B is the leakage coefficient L and B is defined as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðK0=b0Þ=T; Q is the pumping rate, and W(u,r/B) is the leaky well function.
The leaky well function W(u,r/B) can be expressed as W u;r B Z ðN u 1 yexp Ky K ðr=BÞ2 4y dy (3.2)
where y is a dummy variable. Since the right hand side of Eq. (3.2) is an integral form, a numerical approach is required to evaluate the integration. Both the Laguerre quadrature formula and Gaussian quadrature formula (Carnahan et al., 1969) are employed to evaluate the values of leaky well function with the accuracy to the fourth decimal. The Laguerre integration used to approximate an integral function is usually expressed as ðN 0 f ðxÞeKxdx ZX n iZ1 wif ðxiÞ (3.3)
where the wiis weighting factor and xiis
correspond-ing to zero of the nth-order Laguerre polynomials. For a small value of u, the Laguerre quadrature formula cannot give the desired accuracy. Therefore, the Gaussian quadrature formula is employed to evaluate the integration of Eq. (3.2) when u is small.
3.1.2. Dynamic and measurement models for Hantush and Jacob’s model
Parameters T, S, and L, considered as time invariants during the pumping test, constitute the state vector which is to be estimated at each time step The state vector in Eq. (2.2) may be expressed as ^xkðKÞ Z T S L
T
(3.4) After each time step, the renewed state vector, Eq. (3.4), can be substituted into Hantush and Jacob’s model, therefore forming the estimated
drawdown. The estimated measurement from Eq. (2.4) may be combined with Eqs. (3.1) and (3.2) and computed as ^zkZ Q 4pT ðN u 1 yexp Ky K ðr=BÞ2 4y dy (3.5)
3.1.3. Linear approximation equations for Hantush and Jacob’s model
The first-order approximation for the transition matrix for Eq. (2.10) is given as
FKK1z vT vT vT vS vT vL vS vT vS vS vS vL vL vT vL vS vL vL 2 6 6 6 6 6 4 3 7 7 7 7 7 5 (3.6)
Based on Eq. (3.13), Pk(K) can be estimated with
known transition matrix FkK1. Since the state of the
system is described by the parameters, which are mutually independent, thus, FkK1 in Eq. (3.6) is an
identity matrix. Assuming that QkK1is equal to zero,
Eq. (2.9) can be reduced to
PkðKÞ Z PkK1ðCÞ (3.7)
To update the hydraulic parameters in Eq. (2.7), the Kalman gain Kk; estimated by known Hkand the prior
covariance matrix Pk(K), is required. The input zkin
Eq. (2.6) is the drawdown data generated by the Lagrangian polynomial. The measurement matrix Hk
from Eq. (2.8) consists of three elements, which are the partial derivatives of the estimated drawdown ^zkwith respect to T, S, and L, respectively; that is
Hkz v^zvTk v^zvSk v^zvLk
h iT
(3.8) The posterior covariance matrix Pk(C) is then
calculated by inserting the Kalman gain and prior covariance matrix Pk(K). The posterior covariance
matrix is used recursively as the prior covariance at the next time step.
One of the criteria to stop this recursive process is
TOLTZ jTðk C 1Þ K TðkÞj (3.9)
where TOLT is the tolerance criterion for T.
The tolerance criterion TOLs for S and TOLL for
L can be expressed in a similar way. When all
the criteria are met, the recursive process is terminated and the optimal parameters may then be determined.
3.2. Leaky aquifer with storage effect in aquitard Here, the EKF method along with Neuman and Witherspoon’s model (1969) is applied to leaky aquifer considering the effect of aquitard storage Given the initial guesses of parameters, error covariance, and measurement noise, the four best-fit parameters are identified when the convergence criteria are met.
Neuman and Witherspoon presented a closed-form solution for the problem of flow to a well in a confined infinite radial system composed of two confined aquifers that are separated by an aquitard (Neuman and Witherspoon, 1969). Differing from
Hantush and Jacob’ work (1955), Neuman and Witherspoon’s model includes the effect of the aquitard storage on the drawdown of the pumping aquifer. Their mathematical model may be written as s Z Q 2pT ðN 0 1 y½1 K expðKy 2t DÞJ0½wðyÞdy (3.10) where tDZ tDL2=16j2, tDZTt/r 2 S, LZr/B, jZb/L, BZ ffiffiffiffiffiffiffiffiffiffiffiffiffiTb0=K0 p ; bZr ffiffiffiffiS0 p =4BpffiffiffiS; and w2ðyÞZL2y2=16j2 KL2y cot y. Note that Eq. (3.10) is valid for all values of time intervals and the Bessel function J0[w(y)] must be set to zero when w2(y)!0.
The state vector of Eq. (2.2) at each time step is ^xkðKÞ Z T S L j
T
(3.11) Note that those four parameters are also taken as time invariants during the pumping test. After each time step, the renewed state vector ^xkðKÞ from Eq. (3.11) can be substituted into Neuman and Witherspoon’s model, consequently forming the estimated drawdown ^zk: The estimated measure-ment ^zk from Eq. (2.4) is expressed as
^zkZ Q 2pT ðN 0 1 y½1 K expðKy 2t DÞJ0½wðyÞdy (3.12) The measurement ^zk estimated by Eq. (3.12) is used in Eq. (2.6) to obtain the updated state vector. The first-order approximation for the transition
matrix of Eq. (2.10) is given as FkK1z vT vT vT vS vT vL vT vj vS vT vS vS vS vL vS vj vL vT vL vS vL vL vL vj vj vT vj vS vj vL vj vj 2 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 5 (3.13)
The relation between Pk(K) and PkK1(C) can also
be described by Eq. (3.7). Both Hk and Pk(K)
should be calculated in advance by Eqs. (2.8) and (2.9), respectively. The input zk is the measurement
data produced by the Lagrangian polynomial. The measurement matrix Hk is composed of v^zk=vT; v^zk=vS; v^zk=vL; and v^zk=vj; which can be approxi-mated by the finite-difference formula.
The posterior covariance matrix Pk(C) is then
updated by the Kalman gain and the prior covariance matrix PK(K). The updated covariance Pk(C) can
function as the prior covariance in the next time step for the recursive process.
The tolerance criteria of the four-parameter model are the same as those of the three-parameter model, except that the additional criterion, e.g. TOLj, is for the parameter j.
3.3. Evaluation of Hk
InHantush and Jacob’ model (1955), the value of the drawdown depends on the decision variables, T, S, and L, which may be written as
s Z Q
4pTwðT; S; LÞ (3.14)
where the leaky well function w(T,S,L) represents the integral in Eq. (3.2).
The derivatives v^zk=vT; v^zk=vS; v^zk=vL, are, respectively, v^zk vT Z K Q 4pT2w C Q 4pT vw vT (3.15) v^zk vS Z Q 4pT vw vS (3.16) v^zk vLZ Q 4pT vw vL (3.17)
where vw/vT is approximated by forward differencing as
vw vTZ
wðT C DT; S; LÞ KwðT; S; LÞ
DT (3.18)
and the other partial derivatives can also be expressed in a similar manner. The increment shown in the denominator of Eq. (3.18) may be approximated by the parameter value times a factor of 10K3 or less, e.g. DTZ10K3T
In the four-parameter model, the values of draw-down depend on the decision variables, T, S, L, and j, which may be written as
s Z Q
2pTGðT; S; L; jÞ (3.19)
where the function G(T,S,L,j) represents the integral in Eq. (3.10).
In this model, three derivatives v^zk=vT; v^zk=vS; and v^zk=vL are, respectively, given in Eqs. (3.14)–(3.16) and the additional derivative v^zk=vj is
v^zk vjZ Q 2pT vG vj (3.20)
where vG/vT is approximated by forward difference as
vG vT Z
GðT C DT; S; L; jÞ KGðT; S; L; jÞ
DT (3.21)
and the other partial derivatives vG/vS, vG/vL, and vG/vj are also expressed in a similar manner.
4. Two practical problems
4.1. Problems inChander et al. (1981)
Chander et al. (1981) used the Kalman filter to identify the parameters of confined and leaky aquifers. For the leaky aquifer, they combined the Kalman filter with approximate drawdown function without con-sidering the effect of the aquitard storage. Their well function was a truncated form of Eq. (3.2) and may be only suitable for long pumping time or/and small
leakage coefficient. The drawdown equation that they used was s Z Q 4pT 2K0 r L Kw Tt L2S (4.1) This equation may be expanded as
s Z Q 4pT!2 1 C r 2L 2 C 1 4 r 2L 4 C. !ln 1:123r L C r 2L 2 C 3 8 r 2L 4 C. K Q 4pT! K0:5572 Kln Tt L2S C Tt L2S K 1 4 Tt L2S ð4:2Þ The results from applying Chander et al.’s method may have significant errors when the pumping time is short or/and the leakage coefficient is large. In addition, there are two drawbacks in their study. First, laborious works are required to estimate the initial guesses before applying the EKF. The guess values of T and S were determined by using a truncated Theis equation (Chander et al., 1981, Eq. (4)) to the first two drawdown data. The guess value of L was calculated by using the last available observation and applying the approximate steady-state drawdown relationship for the leaky aquifer as used in the Hantush I method (Kruseman and de Ridder, 1970). These strategies may not yield better initial estimates and lose the superiority of on-line estimation of EKF over other parameter estimation methods. Second,Chander et al. (1981)did not specify stopping criteria for their EKF algorithm. It seems to us thatFigs. 1 and 2in their manuscript indicate that the EKF identification process terminated at the time when running out of the observation data. Under this circumstance, the estimated parameters are not warranted to have desired accuracy. The EKF usually takes more than 200 steps to get convergent results, as shown later in this study. In reality, the number of steps needed when using the EKF to determine the hydraulic parameters depends on those guess values mentioned above. The available observed pumping data given in Chander et al.’s study is no more than 12. Thus, the EKF may run out of data very quickly before the hydraulic parameters converge to satisfactory results.
4.2. Problem of small amount of drawdown data Usually, the changes for the estimated para-meters and the elements of error covariance matrix
Fig. 2. The pumping test data ofSridharan et al. (1987)and the estimated drawdown based onNeuman and Witherspoon’s model (1969)for the leaky aquifer with considering the storage effect. Fig. 1. The pumping test data ofCooper (1963)and the estimated drawdown based onHantush and Jacob’s model (1955)for the leaky aquifer without considering the storage effect.
of the state vector in each step of the EKF are fairly small Therefore, the parameter estimation using EKF usually requires the preparing action of a large amount of the drawdown data. However, the number of field measurements is usually limited, say less than 50. Two approaches suggested here to solve this problem are: (1) applying the Lagrangian polynomial to interpolate the field measurement data with smaller time intervals, and (2) recursively using the field measurement data during the step-wise identification process. The first approach is suitable to be used for on-line measurement of the field-pumping test, whereas the second approach is recommended to be applied for an existing drawdown data set.
5. Data analyses and discussion 5.1. Two sets of pumping test data
Two sets of pumping test data are chosen, one for the three-parameter model and the other for the four-parameter model. Three time–drawdown data sets measured from observation wells, as reported in
Cooper (1963) and cited by Lohman (1972, p.31, Table 11), are selected for data analyses. The distances between the pumping well and the obser-vation wells 1, 2, and 3, are respectively, 30.48, 152.4, and 304.8 m. The pumping rate Q is 5450.98 m3/day, the thickness of the aquitard is 30.48 m and total pumping time is 1000 min (16.67 h). The data set for the four-parameter model is taken from Sridharan et al. (1987). The distance between the observation well and pumping well is 29.0 m and the pumping rate Q is 136.26 m3/day. The total pumping time is 1200 min (20.00 h). These two sets of the pumping test data are shown inTables 1 and 2and interpolated using the Lagrangian polynomial.
The time intervals of the pumping test data are normally non-uniform. A small and uniform time interval, say 6 s (i.e. 0.1 min), was chosen when applying the Lagrangian polynomial to interpolate the data. In order to smooth the interpolation results, the time coordinate (t) was transferred to the logarithmic scale (log t) before interpolation and then the second-degree Lagrangian polynomial was employed to interpolate drawdown data for the drawdown versus log t. At the beginning of
the pumping test, the first three drawdown data points were utilized to generate interpolated data between the first and third data points. Once a new drawdown point is measured, the second, third, and the new measured data point (i.e. the fourth data point) were used to generate interpolated data between the third and new data points. Such an interpolation procedure can keep going on in the field until the estimated parameters meet the convergent criteria. For the data sets from Cooper (1963) and Sridharan et al. (1987), the total numbers of the interpolated drawdown data generated by the Lagrangian polynomial are 10,000 and 12,000, respectively. The interpolated time–drawdown data are then ready for EKF identification process.
Table 1
Time–drawdown data obtained from three observation wells (Cooper, 1963, p. 31)
Time (min) Observation well (m)
1 2 3 0.2 0.536 0.003 0.000 0.5 0.838 0.043 0.000 1 1.094 0.137 0.006 2 1.298 0.284 0.043 5 1.609 0.536 0.168 10 1.798 0.713 0.302 20 1.972 0.869 0.445 50 2.109 1.009 0.594 100 2.167 1.067 0.640 200 2.195 1.070 0.643 500 2.198 1.073 0.643 1000 2.198 1.073 0.643
The unit of drawdown is meter.
Table 2
Pumping test data (Sridharan et al., 1987, p. 170)
Time (min) Drawdown (m)
5 0.3 28 0.95 41 1.1 60 1.25 75 1.34 250 1.75 500 1.9 700 1.95 970 1.98 1000 1.99 1200 1.99 QZ136.26 m3 /day; rZ29.00 m.
5.2. Assessment for prediction errors
To assess the accuracy of the estimated parameters or the goodness-of-fit of the estimated drawdown to the observed drawdown, two error criteria, mean error (ME) and standard error of estimate (SEE), are used to calculate the prediction errors between the observed and predicted draw-downs. The ME is defined as
ME Z1 n
Xn iZ1
ei (5.1)
The principle of least-squares assumes that the errors are normally distributed with zero mean and constant variance (McCuen, 1985). When the ME value is equal to or very close to zero, the assumption that errors having zero mean will be satisfied.
The SEE is defined as
SEE Z ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 n Xn iZ1 e2 i s (5.2) where v is the degree of freedom, which is equal to the number of observed data points minus the number of estimate parameters.
5.3. Parameter identification for leaky aquifer without considering the aquitard storage
The initial guesses of the aquifer parameters T, S, and L are given in Table 3. For the initial error covariance matrix, PkK1(C), the diagonal elements
could be assigned as, for example, 55,000, 1.0!10K9 and 1.0!10K6, for T, S, and L, respectively; while and the off-diagonal elements of PkK1(C) are set to
zero. There is no specific rule to decide the initial
Table 3
Initial guess values for EKF and the estimated parameters and prediction errors when using EKF, NLN, and graphical method to analyze Cooper’s data (Cooper, 1963) for leaky aquifer without considering the effect of aquitard storage
Three-parameter model
Case no. Initial guesses for hydraulic parameters Initial error covariance matrix for hydraulic parameters
T S L T S L
1 1000 2.00!10K4 1.00!10K2 55,000 1.00!10K9 1.00!10K3 2 1000 2.00!10K4 1.00!10K1 55,000 1.00!10K9 2.00!10K2 3 1000 2.00!10K4 1.00!10K1 55,000 1.00!10K9 1.00!10K2
Case no. Estimated parameters Prediction errors Step
T S L ME RMSE SEE
EKF on interpolated data
1 1257.9 9.09!10K5 4.82!10K2 K6.53!10K4 1.46!10K2 1.69!10K2 236 2 1311.4 9.29!10K5 2.28!10K1 3.72!10K3 7.46!10K3 8.62!10K3 1508 3 1228.0 1.00!10K4 5.08!10K1 K2.44!10K4 3.42!10K3 4.09!10K3 7603
EKF on recursively used data
1 1239.4 9.78!10K5 4.94!10K2 1.56!10K4 1.15!10K2 1.33!10K2 14,656 NLN 1 1239.1 9.80!10K5 4.94!10K2 K1.10!10K4 1.15!10K2 1.33!10K2 2 1242.1 9.80!10K5 2.52!10K1 4.98!10K4 4.93!10K3 5.71!10K3 3 1215.2 9.70!10K5 5.12!10K1 K1.90!10K4 3.22!10K3 3.72!10K3 Graphical method 1 1234.5 1.00!10K4 5.00!10K2 K2.15!10K4 1.19!10K2 1.38!10K2 2 1234.5 1.00!10K4 2.50!10K1 2.55!10K3 9.10!10K3 1.05!10K2 3 1234.5 1.00!10K4 5.00!10K1 2.33!10K3 4.13!10K3 4.77!10K3
diagonal values. However, since PkK1(C) is equal to
the expected value of squared estimation error E b ð^xkK1ðCÞK xkK1Þð ^xkK1ðCÞK xkK1ÞT c , a reasonable range of initial guess for PkK1(C) can be easily made.
With the computation of EKF method to analyze the measured drawdown data obtained from observation well 1, the values of PkK1(C) are minimized, as are the
estimation errors at each time step. The measurement covariance Rk is set to 10K5 and kept constant
throughout the identification process. Reasonable magnitudes of the measurement error covariance range from 10K4to 10K6. Thus, the median value of Rk is selected in this study. Given the initial guess
values for T, S, and L, along with the initial error covariance matrix, EKF coupled with Hantush and Jacob’s model can identify the hydraulic par-ameters in a step-wise manner. The identification process is terminated when the specified convergence criteria are all satisfied, i.e. TOLT%10K4,
TOLS%10K8, and TOLL%10K6. The parameters,
i.e. transmissivity, storage coefficient, and leakage coefficient, are successfully determined and achieve acceptable accuracy. The estimated three hydraulic parameters are TZ1257.9 m2/day, SZ9.09!10K5, LZ4.82!10K2, and the related SEE value is 1.69!10K2when analyzing the pumping test data in the observation well 1, TZ1311.4 m2/day, SZ 9.29!10K5, LZ2.28!10K2, and the related SEE value is 8.62!10K3when analyzing the pumping test data in the observation well 2, and TZ1228.0 m2/day, SZ100!10K4, LZ5.08!10K2, and the related SEE
value is 4.09!10K3when analyzing the pumping test data in the observation well 3. The comparisons for the results when using EKF, NLN, and the graphical method to analyze drawdown data are listed inTable 3.
5.4. Parameter identification for leaky aquifer with considering the aquitard storage
The initial guesses of the aquifer parameters T, S, and L are listed in Table 4. With the off-diagonal elements set to zero, the diagonal elements of initial error covariance matrix PkK1(C) could have the
values, for example, of 100, 10K10, 10K4and 4.0! 10K7, for T, S, L, and j, respectively. The values of error covariance matrix PkK1(C) can be assigned
similarly to the descriptions in Section 5.3. The measurement covariance Rk is also set to 10K5, and
is kept constant throughout the identification pro-cess. With initial guesses for the hydraulic par-ameters, the application of EKF method to Neuman and Witherspoon’s model can identify the values of those four hydraulic parameters step-by-step.
The four hydraulic parameters estimated by EKF on interpolated data, also shown in Table 4, are TZ22.6 m2/day, SZ1.73!10K4, LZ1.42!10K1 and jZ3.16!10K4; while the related prediction error of the SEE value is 1.37!10K2. On the other hand, the estimated transmissivity, which is the most important hydraulic parameter, by the NLN is 23.3 m2/day and its SEE value is 1.06!10K2.
Table 4
Initial guess values for EKF and estimated parameters and prediction errors when using EKF and NLN to analyze data reported inSridharan et al. (1987)for leaky aquifer with considering the effect of aquitard storage
Method Initial guesses for hydraulic parameters Initial error covariance matrix for hydraulic parameters
T S L j T S L j
50 2.00!10K4 2.00!10K1 5.00!10K4 100 1.00!10K10 1.00!10K4 4.00!10K7
Estimated parameters Prediction errors
T S L j ME RMSE SEE Step
EKF on interpolated data 22.6 1.73!10K4 1.42!10K1 3.16!10K4 1.49!10K3 1.08!10K2 1.37!10K2 839 EKF on recursively used data 23.3 1.64!10K4 1.33!10K1 3.19!10K4 K5.28!10K5 9.04!10K3 1.13!10K2 2001 NLN 23.3 1.65!10K4 1.34!10K1 7.04!10K4 K1.54!10K3 8.50!10K3 1.06!10K2
5.5. Parameter identification while analyzing the pumping data being recursively used
As mentioned above, the EKF method can be applied to analyze the existing pumping data set in a repeated manner In other words, the pumping data can be repeatedly used once the EKF has gone through all the available data. Again, the first data set ofCooper (1963)and the pumping data ofSridharan et al. (1987)
are chosen for parameter estimations using the three-parameter and four-parameter models, respectively. The analyzed results are also demon-strated inTables 3 and 4.
5.6. Discussion
Table 3lists the estimated parameters and predic-tion errors when using EKF, NLN, and graphical method to analyze Cooper’s data (1963) for leaky aquifer without storage effect Compared with NLN and graphical method, EKF on interpolated data can quickly identify the parameters and achieve good accuracy, using only part of the measured drawdown data. The parameters, i.e. transmissivity, storage coefficient, and leakage coefficient, are successfully determined and achieve acceptable accuracy, com-pared with those estimated by NLN as indicated in
Table 3 when analyzing Cooper’s three drawdown data sets (1963). Note that both computer methods, EKF and NLN, give slightly better estimation for the hydraulic parameters judging from the prediction error of the SEE value. The parameters for observation well 1 is quickly determined and the identification process takes less than 236 steps, i.e. 23.6 min of the pumping time to get convergent results. An aquifer test with long pumping time is usually required for leaky aquifers by conventional graphical methods. The total pumping time for the recorded data are 1000 min (16.67 h). However, a long period of pumping time may not be necessary if the EKF is employed. In this case, EKF can save 97.6% of the time required to obtain results from the pumping test. The total numbers of steps for EKF to get convergent values of parameters when analyzing the drawdown data of the observation wells 2 and 3 are significantly larger than that of the observation well 1. This reflects that a larger distance from the observation well to the pumping well requires a longer
observation time. Furthermore, it is found that the estimate leakage coefficient for observation well 1 varies only marginally at the first few observations, but changes dramatically later on during the identi-fication process. This may be attributed to the fact that there is some time lag between the start of pumping and the leakage effect influencing the drawdown.
Table 4 exhibits the estimated parameters and prediction errors when using EKF and NLN to analyze data reported in Sridharan et al. (1987) for the leaky aquifer considering the aquitard storage. The EKF on interpolated data gives results very close to the estimated parameters and prediction errors to those of NLN. There is a total of 839 steps used in the identification process, i.e. 83.9 min of the pumping time. On the other hand, the total pumping time for the observed data are 1200 min (20 h). Clearly, the EKF method can save 93.0% of the pumping test time in this case.
Fig. 1 shows the observed drawdown and the estimated drawdown of those three observation wells based on Hantush and Jacob’s model. This indicates that the estimated drawdowns fit the pumping test data quite well.Fig. 2exhibits that the plot of the observed drawdown and the estimated drawdown based on Neuman and Witherspoon’s model. The figure implies that the parameters estimated by the EKF method are suitable to represent the hydraulic characteristic of the aquifer.
Tables 3 and 4also display the results of estimated parameters and prediction errors while measurement data are recursively used in the EKF identification process using the three-parameter and four-parameter models, respectively. Both tables indicate that the EKF method gives slightly more accurate results, even though many more time steps are needed, due to the measurement data that is recursively used. The state vector defined in Eq. (2.2) is assumed to be noise-free and independent of time, so the state transition matrix, Eq. (2.10), is time-invariant and is an identity matrix. The associated measurement model, Eq. (2.3), assum-ing known error covariance structure, can be evaluated at any time. Thus, the EKF method can be applied for any sampling time, even for analyzing recursively used data. Nevertheless, there are two advantages involved if applying the EKF on the interpolated data with a smaller time interval. First, the measurement of the field-pumping test and the estimation of the hydraulic
parameters can be done with the same on-line measurements in the field. Second, the EKF method usually takes significantly less steps in the identifi-cation process to get convergent results.
Chander et al. (1981)combined the truncated leaky well function with EKF to identify three hydraulic parameters, so the applicable ranges of the pumping time and leakage coefficient L for using their method may be limited.Table 5shows the relative errors of truncated leaky well function versus various values of u and L, indicating that when the leakage coefficient is large, the applicable range of the truncated leaky well function is small if the relative error is restricted to be less than 5%. On the other hand, the relative errors are larger than 5% when uO0.1 and L!0.2. Clearly, Chander et al.’s method may yield significant errors in parameter estimation when the pumping test data do not fall in the appropriate range, as indicated inTable 5.
6. Analyses for model and data measurement uncertainties
In this study, the uncertainties of the measurement data are represented by white noise and temporally correlated noise, and sensitivity analyses for the EKF method are then performed for data with those two types of noise. The model uncertainty in the parameter
estimation is assessed for the case of employing both three-parameter and four-parameter models for ana-lyzing three data sets.
6.1. Measurement data uncertainty
The MATLAB function randn (m, n) with mZ200 and nZ1 is first chosen to generate a realization of white noise (The MathWorks, 1995). The elements in this realization are normally distributed random numbers with zero mean and unit variance. Each element is then multiplied by 4.09!10K3, which is taken from the case of the least prediction error (SEE values) by EKF for observation well 3 in the three-parameter model. Finally, 130 data points are taken from this realization and divided into 10 data sets, so each data set contains 13 elements (i.e. 13 data points). Hantush and Jacob’s model along with the given parameters, say TZ1000 m3/day, SZ10K4, and LZ0.1, was employed to generate 13 synthetic sets of drawdown data. Thus, a set of synthetic drawdown data was generated by simply adding the adjusted noise data to the predicted drawdown data one by one. Accordingly, 10 sets of synthetic drawdown data with uncorrelated noise were obtained.
The original realization of white noise is employed to generate temporally correlated noise. The MATLAB
Table 5
The relative error (%) of truncated leaky well function for various values of u and leakage coefficient L
u L 0.01 0.02 0.05 0.1 0.2 0.25 0.5 1 0.00001 0.08 – – – – – – – 0.00005 0.22 0.22 – – – – – – 0.0001 0.24 0.26 93.57 – – – – – 0.0002 0.26 0.27 0.58 – – – – – 0.0005 0.29 0.29 0.33 27.01 – – – – 0.001 0.33 0.33 0.35 0.16 – – – – 0.002 0.39 0.38 0.39 0.43 37.40 – – – 0.005 0.52 0.52 0.50 0.50 0.51 1.39 – – 0.01 0.74 0.73 0.70 0.66 0.66 0.68 – – 0.02 1.18 1.17 1.12 1.03 0.90 0.87 1.81 – 0.05 2.80 2.78 2.70 2.50 2.09 1.90 1.63 – 0.1 6.43 6.40 6.27 5.94 5.14 4.73 3.36 7.09 0.2 17.18 17.14 16.91 16.33 14.80 13.94 10.27 12.25 0.5 82.79 82.69 82.13 80.62 76.37 73.80 61.09 49.21 1 – – – – – – – –
function Hamming(n) with nZ5 is used to produce five coefficients of a Hamming window (The MathWorks, 1995). The MATLAB function conv(a, b) is applied to convolve vectors a and b (The MathWorks, 1995), where vector a represents the original realization and vector b represents the coefficients of the Hamming window. Algebraically, convolution is an operation which multiplies two polynomials with coefficients containing the elements of a and b. The product of the convolution is adjusted by a factor of 2.88!10K3, which is calculated based on the data variance and the SEE value for the result of observation well 3. The result after the adjustment represents a new realization with temporal correlated noise. Thus, 10 sets of synthetic drawdown data with temporally correlated noise can be formed in a similar manner.
Table 6lists the results of data uncertainty analysis for the case with correlated measurement noise.
The estimated parameter T ranges from 1016.93 to 1042 m3/day, S ranges from 9.22!10K5 to 9.54! 10K5, and L ranges from 8.88!10K2to 9.60!10K2. The SEE values of all cases are on the same order of magnitude.Table 7shows similar results for data with the white noise. Thus, the effect of data with either white noise or temporal noise is negligible in the identification procedure.
6.2. Model uncertainty
The third type of uncertainty listed in Eisenberg et al. (1989) is the conceptual model uncertainty regarding to the geometrical configuration, major features, and boundary conditions. Generally speak-ing, field hydrogeologic information is never known in sufficient detail. Also, the development of a mathematical model usually depends on some
Table 6
The estimated parameters and prediction errors for drawdown data with correlated noises
Case Estimated parameters Prediction errors
T S L ME SEE STEP 1 1033.97 9.26!10K5 0.092 2.05!10K3 1.02!10K2 673 2 1019.79 9.49!10K5 0.095 4.95!10K3 7.62!10K3 1655 3 1023.91 9.39!10K5 0.094 5.39!10K3 8.83!10K3 911 4 1041.93 9.29!10K5 0.089 6.05!10K3 1.20!10K2 300 5 1021.01 9.54!10K5 0.095 3.89!10K3 6.56!10K3 352 6 1042.00 9.31!10K5 0.090 2.72!10K3 1.02!10K2 312 7 1029.79 9.29!10K5 0.092 6.06!10K3 1.02!10K2 346 8 1016.93 9.29!10K5 0.096 7.35!10K3 1.26!10K2 918 9 1020.55 9.47!10K5 0.095 5.61!10K3 8.14!10K3 1173 10 1039.74 9.22!10K5 0.090 5.22!10K3 1.09!10K2 697 Table 7
The estimated parameters and prediction errors for drawdown data with white noises
Case Estimated parameters Prediction errors
T S L ME SEE STEP 1 1033.97 9.24!10K5 0.092 5.71!10K3 1.06!10K2 1033 2 1019.79 9.49!10K5 0.094 5.89!10K3 8.29!10K3 316 3 1023.91 9.15!10K5 0.090 8.60!10K3 1.31!10K2 400 4 1041.93 9.57!10K5 0.096 1.75!10K3 6.30!10K3 496 5 1021.01 9.44!10K5 0.095 2.39!10K3 8.68!10K3 1154 6 1042.00 9.62!10K5 0.096 1.76!10K3 5.01!10K3 1587 7 1029.79 9.31!10K5 0.092 5.53!10K3 9.71!10K3 1293 8 1016.93 9.42!10K5 0.096 2.60!10K3 9.50!10K3 592 9 1020.55 9.65!10K5 0.097 3.42!10K3 5.73!10K3 1132 10 1039.74 9.49!10K5 0.095 5.65!10K3 8.03!10K3 1001
assumptions and/or simplifications. Thus, the selec-tion of a model for describing a target aquifer system is always subject to some degree of uncertainty.
Hantush and Jacob (1955) deals with a confined aquifer underlain by an impermeable bed and overlain by an unconfined aquifer; and an aquitard with a negligible storage exists between those two aquifers. In contrast,Neuman and Witherspoon’ model (1969)
is composed of two confined aquifers that are separated by an aquitard. Both of these models assume that the hydraulic head of the top aquifer maintains constant during the pumping, although Neuman and Witherspoon’s model considers the effect of the aquitard storage on the drawdown of the pumped aquifer.
Three different data sets obtained from pumping tests in leaky aquifers are chosen to assess the difference due to model uncertainty in parameter estimation. The first two data sets, also mentioned above, are taken from Cooper (1963) and Sridharan et al. (1987), and the third data set is taken fromBatu (1998, p. 265). The estimated parameters based on the three-parameter model are listed inTable 3, and those based on the four-parameter model were listed inTable 4. Accordingly, the storage coefficient of the aquitard, S0, can be calculated with known values of j for the four-parameter model; and the hydraulic conductivity of the aquitard K0can be estimated based on the known values of b0and L for both models, as listed inTable 10. ForBatu’ pumping test data (1998), the confined aquifer is pumped at a uniform flow rate of 625 m3/day. The observation well fully penetrating the aquifer completely is located 105 m away from the pumped well, and its drawdown data are given inTable 8. The thicknesses of the aquifer and aquitard are estimated to be 80 and 28 m, respectively. Walton’s type-curve method was used to determine the hydraulic conduc-tivity values of the aquifer and aquitard. Two different match points were selected, point A was arbitrarily chosen and point B was a point on the overlapping part of the curve, as described in Batu (1998). The estimated values of K, K0, and S, by Walton’s type-curve method and the EKF method using the three-parameter model are shown inTable 9. The prediction errors of SEE and ME clearly indicate that the EKF has the better accuracy for the estimated hydraulic parameters than the graphical approach of Walton’s type-curve method.
Table 8
Pumping test data (Batu, 1998, p. 265)
Time (min) Drawdown (m)
0.00 0.000 48.96 0.055 61.63 0.062 75.60 0.066 101.09 0.073 145.44 0.090 192.96 0.098 270.72 0.108 361.44 0.116 449.28 0.122 576.00 0.135 652.32 0.136 767.52 0.137 881.28 0.138 QZ625 m3 /day, rZ105 m. Table 9
The estimated parameters and prediction errors using EKF to analyze data reported inBatu (1998, p. 265)
Method Three-parameter model
Estimated parameters Prediction errors
T S L ME RMSE SEE
Walton’s type—curve A 875 2.7!10K3 0.30 1.26!10K2 1.38!10K2 1.55!10K2 Walton’s type—curve B 899 2.8!10K3 0.30 9.50!10K2 1.06!10K2 1.19!10K2 EKF 907 3.5!10K3 0.34 K2.68!10K3 4.41!10K3 4.97!10K3 Method Four-parameter model
Estimated parameters Prediction errors
T S L j ME RMSE SEE
Again,Batu’s data (1998)is analyzed by the EKF method for the four-parameter model. The estimated values of K, K0, S, and S0 are also reported in
Table 10. The four-parameter model gives slightly better fit than the three-parameter model for Sridharan’s and Batu’s pumping data sets, judging from the prediction errors shown inTable 10. From the viewpoint of curve-fitting, the four-parameter model is more flexible to fit the measurement data than the parameter model. However, the three-parameter model yields fewer prediction errors than the four-parameter model in the case of analyzing Cooper’s data. This discrepancy may be caused by the fact that the top aquifers of the leaky aquifer system for Hantush and Jacob’s model (1955) and Neuman and Witherspoon’s model (1969) are different. The top aquifer in Hantush and Jacob’s model is an unconfined aquifer, whereas in Neuman and Witherspoon’s model it is a confined one. Different geological settings may produce different response of the groundwater flow system while aquifers are under pumping. For the four-parameter model, each estimated value of S0 is at least three orders of magnitude less than the value of S in the same aquifer system for those three data sets as shown in the table. The SEE values for both the three-parameter and four-parameter models are on the same order of magnitude suggesting that the effect of S0 is so small as to be negligible. Clearly, the three-parameter model is a good choice for representing a leaky aquifer system if S0/S!10K3.
7. Concluding remarks
A method using EKF and the Lagrangian poly-nomial is proposed to identify the hydraulic parameters in leaky aquifer systems. Hantush and Jacob’s model combined with EKF can optimally determine the parameters for leaky aquifers without considering the aquitard storage step-by-step in the identification process. Neuman and Witherspoon’s model can also be employed in a similar manner in leaky aquifer with considering the aquitard storage. The drawdown data with non-uniform time intervals from the pumping test are interpolated by the Lagrangian polynomial to obtain the data with a smaller time interval. This interpolation approach facilitates the implementation of the EKF to estimate the hydraulic parameters on-line in a field-pumping test.
The proposed approach can quickly identify the parameters, using only a few observed drawdown data, and the obtained parameters are shown to achieve good accuracy. At present, pumping tests are usually performed in the field with a system having pressure transducers installed in the obser-vation wells to measure the water level and a data logger to store the measured data transmitted from the pressure transducers. Such a system may be linked to a computer in which the EKF coupled with the aquifer model is implemented and executed simultaneously. Accordingly, the hydraulic parameters can be deter-mined on-line in the field. Once stable estimates of the hydraulic parameters have been reached, the pumping test may be terminated.
Table 10
Estimated hydraulic conductivity and storage coefficient for aquifers and aquitards
Model Estimated parameters Prediction errors
T K K0 S S0 ME SEE
(Cooper’s data (1963, p. 54)a
Three-parameter 1257.9 N/A 9.59!10K2 9.09!10K5 K6.53!10K4 1.99!10K2
Four-parameter 1242.7 N/A 7.89!10K2 1.16!10K4 2.95!10K9 K1.81!10K3 5.79!10K2
Sridharan’s data (1987, p.170)
Three-parameter 21.60 N/A N/A 1.69!10K4 K1.95!10K2 2.73!10K2
Four-parameter 22.62 N/A N/A 1.73!10K4 2.76!10K10 1.49!10K3 1.37!10K2 Date reported inBatu (1998, p. 265)b
Three-parameter 907 11.34 0.26 3.5!10K3 K2.68!10K3 4.97!10K3
Four-parameter 893 11.16 0.26 3.43!10K3 1.31!10K6 K1.51!10K3 4.61!10K3 N/A represents that the data are not available because the aquifer thickness (b) or/and the aquitard thickness (b0) is unknown.
a b0
Z30.480 m.
b
The analyses on Chander et al.’s truncated leaky well function show that their approach of using EKF to estimate the hydraulic parameters is only suitable for cases where the pumping time is large or/and the leakage coefficient is small. The identification process of EKF has been shown to be able to reflect the physical nature of the leaky aquifer. The slow response of the leakage coefficient may be attributed to the fact that there is some time lag between the start of pumping and the leakage effects influencing the drawdown.
Finally, the sensitivity test for the measurement noise indicates that the effect of both the white noise and the temporal correlated noise in measured draw-down is insignificant when employing the EKF to identify the hydrogeologic parameters of the leaky aquifers. The choice between Hantush and Jacob’s model and Neuman and Witherspoon’s model for representing the leaky aquifer system depends on the hydrogeological condition of the system indicated in the analyses of the model uncertainty. However, Hantush and Jacob’s model is suggested for use if the ratio of the aquitard storage to the aquifer storage is less than 10K3.
Acknowledgements
This study was partly supported by the Taiwan National Science Council under grant NSC90-2621-Z-009-003. The authors would like to thank the two anonymous reviewers for their valuable and con-structive comments.
References
Batu, V., 1998. Aquifer Hydraulics. Wiley, New York.
Beck, M.E., 1987. Water quality modeling: a review of the analysis of uncertainty. Water Resour. Res. 23 (8), 1393–1442. Bierkens, M.F.P., 1998. Modeling water table fluctuations by means
of a stochastic differential equation. Water Resour. Res. 34 (10), 2485–2499.
Bierkens, M.F.P., Knotters, M., Hoogland, T., 2001. Space–time modeling of water depth using a regionalized time series model and the Kalman filter. Water Resour. Res. 37 (5), 1277–1290. Cahill, A.T., Ungaro, F., Parlange, M.B., Mata, M., Nielsen, D.R.,
1999. Combined spatial and Kalman filter estimation of optimal soil hydraulic properties. Water Resour. Res. 35 (4), 1079–1088.
Carnahan, B., Luther, H.A., Wilkes, J.O., 1969. Applied Numerical Methods. Wiley, New York.
Chander, S., Kapoor, P.N., Goyal, S.K., 1981. Aquifer parameter estimation using Kalman filters. J. Irrig. Drain. Div., Am. Soc. Civ. Eng. 107 (IR1), 25–33.
Christiaens, K., Feyen, J., 2001. Analysis of uncertainties associated with different methods to determine soil hydraulic properties and their propagation in the distributed hydrological MIKE SHE model. J. Hydrol. 246, 63–81.
Cooper Jr., H.H., 1963. Type curves for nonsteady redial flow in an infinite leaky artesian aquifer, in: Bentall, Ray, complier, Short-cuts and special problems in aquifer tests. US Geol. Surv. Water-Supply Pap. 1545-C 1963;.
Dawson, K.J., Istok, J.D., 1991. Aquifer Testing: Design and Analysis of Pumping and Slug Tests. Lewis Publishers, Chelsea, MI.
Eigbe, U., Beck, M.B., Wheater, H.S., Hirano, F., 1998. Kalman filtering in groundwater flow modeling: problems and prospects. Stochastic Hydrol. Hydraul. 12, 15–32.
Eisenberg, N.A., Rickertsen, L.D., Voss, C., 1989. Performance assessment, site characterization, and sensitivity and uncertainty methods: their necessary association for licensing, in: Buxton, B.E. (Ed.), Proceedings of the Conference on Geostatistical, Sensitivity, and Uncertainty Methods for Ground-Water Flow and Radionuclide Transport Modeling, CONF-870971. Battelle Press, Columbus, OH.
Freeze, R.A., James, B., Massmann, J., Sperling, T., Smith, L., 1992. Hydrogeologic decision analysis. 4. The concept of data worth and its use in the development of site investigation strategies. Ground Water 30 (4), 574–588.
Gerald, C.F., Wheatley, P.O., 1994. Applied Numerical Analysis, fifth ed. Addison-Wesley, Reading, MA.
Grewal, M.S., Andrews, A.P., 1993. Kalman Filtering: Theory and Practice. Prentice-Hall, New Jersey.
Hantush, M.S., 1960. Modification of the theory of leaky aquifers. J. Geophys. Res. 65 (11), 3713–3724.
Hantush, M.S., Jacob, C.E., 1955. Non-steady radial flow in an infinite leaky aquifer. Trans. Am. Geophys. Union 36, 95–100. Jacob, C.E., 1946. Radial flow in a leaky artesian aquifer. Trans.
Am. Geophys. Union 27, 198–208.
Katul, G.G., Wendroth, O., Parlange, M.B., Puente, C.E., Folegatti, M.V., Nielsen, D.R., 1993. Estimation of in situ hydraulic conductivity function from nonlinear filtering theory. Water Resour. Res. 29 (4), 1063–1070.
Kruseman, G.P., de Ridder, N.A., 1970. Analysis and evaluation of the pumping test data, Bulletin 11. Institute for Land Reclamation and Improvement, Wageningen, The Netherlands pp. 80–94.
Lee, Y.P., Chin, C.V., Yen, P.H., 2000. Application of Kalman Filter for Real-time Identifying Aquifer Anisotropy, in: Proceedings of the 11th Hydraulic Engineering Conference, Taipei, Taiwan pp. D81–D86.
Lohman, S.W., 1972. Ground-water hydraulics. US Geol. Surv. Prof. Pap. 1972, 708.
McCuen, R.H., 1985. Statistical Methods for Engineers, Prentice Hall, Englewood Cliffs, New Jersey.
Neuman, S.P., Witherspoon, P.A., 1969. Theory of flow in a confined two aquifer system. Water Resour. Res. 5 (4), 803–816. Sridharan, K., Ramaswamy, R., Lakshmana Rao, N.S., 1987.
Identification of parameters in semiconfined aquifers. J. Hydrol. 93, 163–173.
The MathWorks, 1995. The Student Edition of MATLAB, Version 4, User’s Guide. Prentice-Hall, Englewood Cliffs, NJ. Van Geer, F.C., Van Der Kloet, P., 1985. Two algorithms for
parameter estimation in groundwater flow problems. J. Hydrol. 77, 361–378.
Van Geer, F.C., te Stroet, C.B.M., 1990. A Kalman filter approach to the quantification of the reliability of a groundwater model, in: Kovar, K. (Ed.), Calibration and reliability in groundwater modeling ModelCare 90, IAHS Publication No. 195, pp. 467–476. Van Geer, F.C., te Stroet, C.B.M., Bierkens, M.F.P., 1990. Groundwater modeling in relation to the system’s response
time using Kalman filtering, in: Gambolati, G., Rinald, A., Brebbia, C.A., Gray, W.G., Pinder, G.F. (Eds.), Proceedings of the Eighth International Conference on Computational Methods in Water Resources and Surface Hydrology, Part A. Compu-tational Mechanics Publication, Venice, Italy, pp. 23–30. Wheater, H.S., Tompkins, J.A., can Leeuwen, M., Butler, A.P.,
2000. Uncertainty in groundwater flow and transport model-ing—a stochastic analysis of well-protection zone. Hydrol. Process. 14, 2019–2029.
Yeh, H.D., 1987. Theis’ solution by nonlinear least-squares and finite-difference Newton’s method. Ground Water 25 (6), 710–715.
Yeh, H.D., Han, H.Y., 1989. Numerical identification of parameters in leaky aquifers. Ground Water 27 (5), 655–663.
Zheng, C., Bennett, G.D., 2002. Applied Contaminant Transport Modeling, second ed. Wiley/Interscience, New York.
