- Email: [email protected]
Parameter identification in radial flow
Parameter identification in radial flow R. T. Hudspeth, R. B. Guenther*, K. L. Roley, and W. G. McDougal
Department oj" Civil Engineering. Oregon State University. Corvallis. Oregon 97331. USA The classical inverse parameter identification problem is solved using an indirect method. This method is based on the minimization of an objective function or error criterion consisting of three parts: 1) least-squares error of head residuals; 2) prior information of flow parameters: and 31 regularization. An adjoint equation is incorporated into the method to eliminate the need to differentiate the heads with respect to the parameters being identified, increasing the stability of the algorithm. This method is applied to a radial coordinate system used for describing an anisotropic confined aquifer. A variable implicit finite difference method is used for predicting hydraulic heads given a continuous set of flow parameters. The method is tested using synthetic data and is found to give estimates which are both stable and unique.
1. I N T R O D U C T I O N
using only values of head according to
During the past 25 years, a number of algorithms have been developed to solve the classical inverse problem for groundwater flow. Yeh 3"*and Willis and Yeh 32 reviewed these algorithms in detail. Neuman e° has proposed classifying parameter identification algorithms as either direct or indirect. The direct algorithm is a non-iterative technique where the parameters are unknowns and the head values are knowns. Since piezometric head measurements are usually available only at a few specified locations in the field, some means of interpolation must be used to estimate heads at other spatial locations resulting in interpolation errors. The results obtained by the direct algorithm also seem to be correlated to the grid size used {Willis and Yeh32). The indirect algorithm uses an iterative technique to minimize an objective function. The least squares error between the measured and estimated heads provides the basis for many objective functions (Neumana'; Cooley~; Sun and Yeh-'r; Carrera and Neuman3"*'s). Additional penalty terms mav be added to the least squares error resulting in a multiple objective function. The least squares criterion is capable of filtering out some of the measurement errors inherent in field data. However, this term alone is not enough to guarantee a unique set of parameters (Carrera and Neuman3). A unique set may be computed by the addition of a prior data penalty term to the objective function. The inverse parameter identification problem may be conceptualized as follows: Given the fvnctional relationship
h(r, t ) = For [p(r)]
(1)
where h(r. t) is the hydraulic head; r and t are space and time, respectively; and plr) are spatially varying parameters that control the groundwater flow: identify p(r) * Department of Mathematics, Oregon State University, Corvallis, Oregon 97331, USA C 1991 Elsevier Science Publishers Ltd.
240
Adv. Water Resources, 1991, Vol. 14, No. 5
p(r) = Inv [h(r, t)]
(2)
Carrera and Neuman'* identify (1) as a "forward relationship' and (2) as an "inverse relationship.' This inverse problem is well-posed in the sense of Hadamard 15 if: (1) For every value of head data h(r, t), there is one and only one solution for the parameters p(r), and (2) The parameters p(r) depend continuously on the head data h(r, t) (John 17, Guenther and Leer'*). Condition (1) insures that the parameters are identifiable and unique. Identifiability is not achieved in the forward relationship if different sets of flow parameters result in the same head distribution. Uniqueness is not obtained in the inverse relationship when different head distributions produce the same set of flow parameters. This distinction is important since an indirect parameter identification algorithm may produce an identifiable solution that is non-unique (Carrera and Neuman'*). Condition (2) insures that the values of the parameters are stable. The parameter values are unstable if small errors in the head data result in large errors in the parameter values. Two types of errors are known to cause instability: 1) modeling and interpolation errors, and 2) measurement errors (Neuman and Yakowitzee). Analyses of the sensitivity coefficients indicate that errors in the head data contribute the most to instability. The sensitivity coefficients are defined as the change in head with respect to the change in the parameters being identified. Values of sensitivity coefficients that are nearly zero indicate for the inverse relationship that a small change in measured head (i.e. measurement error) may result in a large change in the parameter values (Yakowitz and Duckstein33).
1.1 htentifiabilit.v Sufficient values of measured head data are assumed to be available so that a measurable head gradient exists between any two spatial locations within the flow region during a finite period of time for the duration of the test (Carrera and Neuman'*). It is not always possible to satisfy
Parameter identf/ication in radial flow . Q~(~)
R. T. Hudspeth et al.
.*P~ezomefer
~.-:..:.~-~.: :."..:~.,.q ~ : ~ , T X ~ Z ,
j [:.,.....c..;:.,_..,'.'..'.'.,.
/
".
r=
"
:h =
:t =
" Sir)=
• T(r)=~(3a-e)
E
[\
/
h° {r.lLll
'/
h(r.l~.l) I ~".~
.~," ".7
L(r. ll.l) "\\~'\~x
Flow
Oomo,n
'\ " "
"\
i
The parabolic diffusion equation governing radial flow through porous media shown in Fig. I is well-known (McWhorter and Sunada ~9) and is given in dimensionless form by ~ q h ( r , t)) =
~N'.
F~,. 1.
\ \\
S(r)
~h(r, t)
?t
1 ~ ( (:h(r, t)) Pr ,, rT(r) ~ = 0; ('r
r
\"
r" < r < R . t > ~ O
Dqfinitio~z sketch
this assumption on the boundary of the flow region: and thus the identifiability of the parameters cannot always be guaranteed. Dimensionality may be reduced by estimating the "true" parameter distribution with functions having finite number of coefficients (Carrera and Neuman4). 1.2 Uniqueness
(4)
in which ~'(. ) is a linear operator, S(r) is storativity, T(r) is the transmissivity, r" is the well radius, and R is the radial distance to the boundary of the flow domain. The initial and boundary conditions may be scaled by riO(r) /~h(t) h°(r) = Ft~ " H a ( t ) - f'f (~'(0. Qh(t)" Q'(t) = ~ , Qb(t) = ~ ;.
)~ T~
(5a-e)
The dimensionless initial and boundary conditions are
Unique parameter values result from the indirect method when prior data are included in the error criterion (Neuman2°). Neuman et a[. 23, Guenther, eta[. ~3, Carrera and Neuman 4's and Gerlach and Guenther ~ demonstrate the use of prior data. Prior data are parameter values that have been "measured" or estimated at specific locations within the flow domain. Prior data may be obtained from eflher field or lab tests or by comparison with similar aquifers. It must be assumed that all measured parameter values contain some noise, the magnitude of which may or may not be known. 1.3 Stability Stability may be obtained by several methods. The adjoint equation (Chavent 6) may be used to eliminate all sensitivity coefficients from the gradient equations used to minimize the objective function. The relationship between instability and the sensitivity coefficients is demonstrated by Yakowitz and Duckstein 33. Some algorithms that use the adjoint equation do not eliminate the sensitivity coefficients from the gradient equations used to minimize the objective function (Sun and Yeh2:). An additional penalty term composed of the head residuals at the final value of time may be added to the error criterion. The addition of this penalty term improves stability by taking advantage of more reliable values of measured final heads due to repeated measurements and the reduction in the rate at which the head values are falling. Finally, Kravaris and Seinfeld ~s and Gerlach and Guenther ~ add a regularization term to the objective function to improve stability.
h(r, 0) = H°(r), lira
rT(r)
~h(r, t) ?r
rT(r) -?h(r, t) ~r
-
r"'< r < R
(6a)
-Q'(O
(6b)
'
t > 0
2~
,;. 1 27z h(r, t) = ~ [--,;.Hb(t) -- Qh(t)]; r = R, t > O
(6c)
where tfh(t) is the prescribed head at the radial boundary; Qh(t) is the prescribed flux across the radial boundary; and ). is a leakage coefficient at the radial boundary. Note that the leaky, head-driven boundary condition given by (6c) may be used to prescribe Dirichlet boundary conditions when ). = :,-. and Neumann boundary conditions when ). = 0. Numerical solutions may be improved if the arbitrary time-dependent pumping function Q'(t) is defined by N
Q"(t) = y" [Q~ + Qii(t)].U,it. r, ,, r,,)
(7)
n=I
where N is the total number of finite time intervals during the total test duration; Q~, is a constant pumping rate during time interval n; Q~(t) is a variable pumping rate during time interval n: and the Heaviside Step function U,(t, ~,,_ ~, L,) is a pathological numerical switch defined as (Cochran et a/.7), U,(t, v,_ 1, ~,) = U(t - z,_ 1 ) - U(t - z,)
(8)
that turns the well pumping function Q'(t) on at t = v,_ and off at t = % during the time interval n. 2. la Hybrid analytical-numerical solution
2. B O U N D A R Y
VALUE
PROBLEM
(BVP)
Dimensional variables denoted by (?) are scaled by dimensional scaling parameters denoted by (')s to obtain dimensionless variables according to
The strategy used in the hybrid analytical-numerical approach is based on linear superposition (Stakgold26). The solution to the dimensionless BVP (4 - 6) is assumed to be a linear combination of the following linearly independent functions:
Adv. Water Resources, 1991. Vol. 14, No. 5
241
Parclme[cr idcnHticatiim tn rcldial jlon. R. T. HuJ~pcth et al. hit. ti = h ] [ I C ] :
[BC~]: [BC:]: [I-]', _
= h,,:[IC]:
0: 0: 0', + hi',0: [BC,]:
0: 0',
where [IC] are initial conditions: [-BCi] are boundary data at the ith boundary': and [ j ] is the forcing function. The linear superposition principle suggests looking for a linear combination of linearly independent functions which satisfy only specific elements of the inhomogeneous system while contributing nothing to the other inhomogeneous elements (i.e. they' satisfy homogeneous conditions elsewhere). In principle, a linear superposition approach would seek the linear combination (9b): however, in order to take adwmtage of the well-known Theis integral (Theis es, Guenther and Lee~'~), the ideal linear combination (9b) will be replaced by' (I0)
where h°(r, r) satisfies the initial conditions: and ,4r, t) and r(r, t) are analytical and nurnerical approximations to the b o u n d a r y data. Substituting (10)into (4) gives 5/'(h) = 5/'(h °) + 5/'(u) + 5/'(r) = 0
(l I)
The initial data are satisfied exactly by the solution to 5/'(h°): while the b o u n d a r y data are satisfied by the solutions to 5/'(u) and L/'(v). The following dimensionless BVP is defined for t/~(r, t): L/'(IC) = S(r) Fh°(r, t)
I C IrT(r) (h°(r, = t) 1
(:t
=0:
r" < r <
r 7r
(12a)
R,t>0
with inhomogeneous initial data h°(r. 0)= H°b'):
r'"
R
(12b)
and homogeneous boundary data ~:h°(r,
t
0:
c~r
r= R:t >0
(12c)
The following dimensionless BVP with constant coefficients is defined for u(r, t):
~,O(u)=So&4r.
t)
(':[
T° F (Cu(,', t)) -
=
I"
r
.
(T
( F
=0:
,
r'" < r < R, t > 0
(13a)
with homogeneous initial data u(r, 0 ) = 0,
r " ~ < r ~< R:
(13b)
with inhomogeneot, s b o u n d a r y data at the well lim r
r' ~ 0
rT o ?u(r, t) dr
-
-Q"(t) 2=
t >
O:
(13c)
and with a homogeneous far field condition u(r,t)=0.,
r> R,t >0
(13d)
The parameters S o and T o are values of these parameters at the well r = r" --+ 0. The dimensionless solution to BVP (13) is the well-known Theis integral given by
242
Adv. Water Resources. /991. Vol. 14. No. 5
t't
o O- " ( t ' )
S(r)
Ft'(r, t) ~t
l ((rT( v r- ~r
(14)
(r - t)
&'(r, t)) r) ~ : T f - - ,
. =
l(r,
t)
r" < r < R , t > 0
(15a)
t)] -~:u(r" - - - t) + I .F ("rT(r)(U(r. . . .
(15b)
with forcing function J(r,t)=-S(r)
~t
r ~r
~'r /"
with h o m o g e n e o u s initial condition dr, 0 ) = 0
r" ~ r ~ < R
(15c)
with inhomogeneous b o u n d a r y data at the well lira r=r"
rT(r)
FRO', t)
-
~h°(r, t) rT(r) - -
t > 0
(15d)
~0
and with inhomogeneous far field b o u n d a r y data ~:v{r,t)
rW(r)
Fr
}.
1
2= v(r,t) = ~ - rT(r)
[
}.Hh(t) -- Qh(t)]
~:u(r, t) " r = R. t > 0 ~r
(15e)
The solution for h(r, t) is now defined by the following linear system of dimensionless BVPs: h(r, t) h°(r, O)
{[IC]: [BC,]: EGG] [0]I
~r
2
The dimensionless BVP for r(r, t) becomes
(9a,b)
h(r, t) = ttr~(r, t) + u(r, t) + dr, t)
i
ulr, t) - 4]4-7%
+ hj.',O: 0: 0: [./]i
It_,',(): 0: [BC:]:0;
-
Sot
+
u(r, t)
',[o]: [BC, ]: [0]: [0]',
+
',E/c]: [o]: [o]. [o]', r(r, t)
',[o]: [BC,]; [~C_,]; [.I-]I
The Theis integral was c o m p u t e d using both infinite series and G a u s s - L a g u e r r e quadratt, re (Roley2a). The forcing function f(r, t) in { 15b) must be calculated from u(r, t). The value of f(r, t) represents the residual to (ll) when the distributed parameters T(r) and S(r) are not constants. If T{r) and S(r) are constant, then f(r, t) =- O. The hybrid analytical-numerical method offers several advantages over the totally numerical method. (1) Lower sensitivity to At. During the initial interval of the p u m p test, the analytical solution by the Theis integral (14l is independent of At and accounts for most of the drawdown. Therefore, regardless of the p a r a m e t e r distribution, the initial solution to (I 1) is nearly analytical with r(r, t) --- 0. This is a particular advantage because d r a w d o w n is most rapid at the initiation of pumping, and numerical errors tend to be greatest during this interval. (2) Smaller numerical errors. The numerical solution r(r, t) is only' required to correct the Theis integral (141 when the parameters are non-constant or distributed. Sensitivity' to numerical methods is thereby reduced. (3) Separation of variables. The linear superposition given by {16) provides a convenient way to separate the effects on the Theis integral of an aquifer having leaky b o u n d a r y conditions, distributed parameters, etc.
Parameter identification in radial flow. R. 7". Hudspeth et al. A disadvantage of the hybrid analytical-numerical method is an increase in the computational time compared with a fully numerical method. For example, a solution using 75 spatial nodes and 50 time steps on a personal computer required 140% more computer time.
2. lb Numerical method The radial derivative ~h(r, t)/& in (6b) approaches infinity at the well as a result of a singularity at r = 0. To approximate this derivative accurately near the well boundary in discrete form, a very small Ar is necessary. A coordinate transformation that was introduced by Ehlig and Halepaska 9 and used by Ruston and Redshaw 25 maps a variable Ar into a constant At) by
for transient data used by Carrera and Neuman ~'. Both of these objective functions include weighted least squares error and prior data terms. The objective function 118) differs from most other objective functions by the addition of the following terms: (1) A weighted least squares error in the measured and predicted heads, Jflp), at the final value of time, r x, which optimizes the improvements in data measurements made at later times. (2) A regularization term, Jg(p), which smooths the parameter estimates if errors are present in the measured data.
The dimensionless BVP (4 - 6) in r is transformed to
In addition to adding these two terms to the objective function, the BVP ( 4 - 6) is scaled in order that each error or parameter estimate in the objective function are of the same order of magnitude. This results in a faster convergence rate for the parameter identification algorithm.
(r"2exp[2p]S(P) g'h(p'
3.1 Least squares error, JnfP)
[)p~ [T(p) ch(p' t)] = OO
(17a)
A weighted least squares error between the "measured' or field heads, h"(r, t), and the predicted heads, hP(r, t), is defined by
{17b)
Jh{P) = ~
with dimensionless initial conditions
h(p,O)= H°(p);
0
P
lira T(pJ ,,=~
T(p)
?h{p, t) ?p
~h(p. t) (p
- Q"(t) - - - " 2r~ '
Vl(r, t)[hP(r, t) - h (r,t)]-rdr;
dt L
and dimensionless boundary conditions
r"
1 ~n<~N t>0
(17c)
) h(p, t)= 1 27z ~ [-- )'Hb[t) - Qb{t)]; p= P,t>0
(17d)
where P = In (R/r'") and ?h(p, t)/@ at p = 0 is now finite. The dimensionless transformed BVP (17) may now be solved very etticiently by numerical methods.
2. lc Numerical solution A numerical solution was computed using a weighted implicit finite-difference method (Carnahan, et al.2). The value of the weighting factor 0 may be chosen such that when 0 = 0, a fully explicit method results; when 0 = 0.5, the Crank-Nicholson method results; and when 0 = 1, a fully implicit method results. Numerical results indicated that 0 = 0.6 provided the best estimates. At p = P, a fictitious node is defined outside the domain using the boundary simulation technique described by Hildebrand ~6. The boundary condition at p = P is then used to eliminate the fictitious node. The same technique is also applied at the well, p = 0 (Roley-"~). The dimensionless BVP (17J is reduced to a set of algebraic equations in which the coefficient matrix is tridiagonal. The Thomas algorithm (Vemuri and Karplus 3t) was used to solve this matrix. 3. O B J E C T I V E F U N C T I O N The objective function consists of four terms defined as (Carrera and Neuman "t)
J(pt = Jh(P) + Ja(P) + Jf(P) + JR(P)
(18)
The objective function (18)is similar to equation (24)
(19)
where V~(r, t) is a weight which represents the confidence in the measured head data h"(r, t). The least squares error is the basis for most indirect parameter identification algorithms. However, this term alone is not enough to insure either stability or uniqueness in the parameter estimates.
3.2 Prior data error. Ja(P) In order to achieve uniqueness, a prior data error is defined by 1 ~
Jd(P) = ~ L -
,
V2t[Se(~t) -- sm "(~t)]-
/=I
+S
1
V3t[Tt{tlt) - TmOhl] 2
(20)
- l = I
where S"(~ t). . . . . S"(~.,t) and T"(qt) . . . . . T"(qr) are measured values of storativity and transmissivity, respectively, measured at specific locations within the confined aquifer (i.e., ~l . . . . . ~.~t and ql . . . . . qr); and Se(~l) . . . . . Se(~,,t) and Te(q~) . . . . . T~(qr) are the estimated values of storativity and transmissivity, respectively, obtained from a parameter identification algorithm. The weights V2t and V3~ represent the confidence in the measured prior data. The addition of prior data may not reduce the least squares error observed in the head values. Because of possible errors in the measured head data and the difficulty in obtaining accurate estimates of the flow parameters from field data, the minimum value of the objective function may even increase as a result of adding prior data (Neuman2~). Previous attempts to compute unique parameter values based on experimental data and/or professional judgement have yielded reasonable parameter estimates. The two remaining terms in the objective function (18) are added to provide stability.
Adv. Water Resources. 1991, Vol. 14, No. 5
243
Parameter ident(/ication in radial flow. R T. Hmtspeth et al. 3.3 Final head error. Jr, p,,
The heads measure at the final value of time L~ may be more reliable than earlier measured values because of repeated calibration checks and improved data recording procedures. Accordingly, these final measured heads are weighted more heavily. If the reliability of measurements at a particular location r~ are uncertain, the weighing coefficient l:~.(r~)for these locations may be set to zero.
descent and Fletcher-Reeeves gradient methods (Bazarra and Shetty~). First. the adjoint BVP is derived using variational calculus to eliminate the sensitivity coefficients from the gradient equations. Then. the steepest descent and Fletcher-Reeves discrete gradient methods are used to minimize the objective function. Two potential difficulties ma? be encountered when calculating the gradient equations. First, these gradient equations include sensitivity coefficients which may lead to instabilities. Second, these sensitivity coefficients must be calculated from a BVP requiring the solution of N t + N, + l differential equations. The adjoint BVP may be used to eIiminate both of these potential difficulties.
3.4 Tikhom~v regularization. J~(p)
4.1 .q~Ooittt B I'P (Backward Heat Equation)
Since instabilities in the values of parameter estimates appear to be closely' related to noise in the measured data, a Tikhonov regularization term is added that smooths the parameter estimates. Although there are several ways this could be done, the method of Tikhonor :~''~° used by, Kravaris and Seinfeld Is and by Gerlach and Guenther ~1 is used and is defined by
The adjoint BVP may be derived by first multiplying the variational of the state equation (4) by a scalar function w(r, t), and then integrating by parts over both the spatial and temporal domains to obtain (Friedman ~°)
A weighted least squares error between the measured and the predicted heads at the final value of time is defined bv I
" R
J /(pl = ~_..... Va(r)Slrl[hPlr. r~.) -- h"(r, r~.)]2rdr
= 2
+
121)
S,.,,.) ~',,'(,', - - t) + -1 ( f-~. )
?t
Jrdr
(22)
where Vy(r)and l/~,(r)are weights. This regularization term imposes a penalty on oscillations in parameter estimates by increasing the value of the objective function (Carrera and Neuman3l. The use of the adjoint equation provides an important contribution to stability of the parameter identification algorithm. The adjoint equation can be used to eliminate the need to differentiate the measured heads with respect to the parameters being identified. Note that all of the error terms, Jh(p), Ja(P), and J/(p), reduce to zero when the predicted and measured values are identical. However, for the regularization term (22), J~(p) is always positive-definite. This is not a problem since the objective function is minimized relative only to an initial value. This point is important when establishing criteria for terminating iterations in a parameter identification algorithm.
3.5 Reduction of dimensionality of parameters The distributed parameters, T(r) and S(r), may be approximated by constructing sequences using simple polynomial basis functions according to
j=l
,,,,.flrI" I~
4. M I N I M I Z A T I O N FUNCTION
The objective function is minimized by both the steepest
244
Adv. Water Resources. 1991. Vol. 14. No. 5
~r
(25a)
with inhomogeneous initial data w(r, r x ) = -l,;,(r)[hQr, r s) - h " ' ! r , rs)]:
r'~
and homogeneous boundary data lim
rT"lr)
~;w(r, t)
-0:
0
(25c)
2re w(r. t): 0 < t < r x, r = R
(25d)
~2r
r=r"~O
and
rTqr)
~w(r, t)
(r
)."
where ): is a leaky parameter that must also be identified. If V dr) =- 0 in (25b), the homogeneous initial condition used by Carrera and Neuman "~5 is obtained. The solution to the adjoint BVP (25) reduces the first variation of tile objective function (18)to
dt ill ~hP("" = t)(~,c[,-, t)
J(.)(p) =
(T
~~r
dt -R Hip(r, ( ~ - t)
w(r,
j
(23a,b)
OF THE OBJECTIVE
I rTqr) ,'w(,', t)]
r'"~r~R,O
m=l
where the estimates Tqr) -+ T(r) and S"(r)---, S(r) as N i, N_, ~ m. Parameter identification is now reduced to identifying a finite number of coefficients, :~j and tim, in (23) such that the objective function (18) is minimized.
r ~r
= V~(r, t)[/f'(r, t ) - h"(r, t)]:
N.
Tqr)= ~. :~j"r Ij 1~., S~'(r)= 2
(24)
where ~'*(.) is the adjoint operator. The inhomogeneous adjoint BVP for w(r, t) is
w
/~Tqr)~2q
N ~
T~',O')rdr
[)S[.)(rt
rdr
r~
.If
+E I-I
-
K
+E I
~[~t[T"(,l,)
-
T"(r/,)]
T~I~(,1,)
(26)
1
l:~(r)Se(r)[hQr, r x) - h " ( r , rx)]es[l(r)rdr
Parameter ident(/ication in radia/ ffou. R. 71 Hudspeth et al. "R
- t
-
-
[drlS~'(rIS"lrwdr [~lr) T~lr)T~',lr)
r"
(T"lr)= (T~"lr)]rdr (r
+
i
,:.~',[Hh(t}-- hPlr,tl]wlr, Hdt
where JcjIp) = ~ J { p ) ~ { ")is the partial derivative of(18) with respect to the parameters :~ and fl in {231 and 2" in 125d). Note that the only head derivatives which are to be computed are those of the predicted head values hP(r, t) which are smooth. Derivatives of the measured head values (the sensitivity coel~cients) have been completely eliminated from the gradient equations. Partial derivatives with respect to the parameters to be identified are required to minimize the objective function {18). In both of the gradient methods {steepest descent and Fletcher-Reeves), a one-dimensional line search procedure was employed for each iteration to determine a local minimum of the objective function for each parameter group along their respective directional vectors. A global optimum is defined as a set of all of the parameters which minimize the objective function with respect to the minimization criteria.
4.2 Steepe,vt ek',s'cent gradient method The steepest descent method is one of the oldest and simplest gradient methods used (Bazarra and Shetty~). The optimum step size is a function of the parameter being identified. Selecting a single step size for all the parameters results in less computer time per iteration. Unfortunately, this often makes convergence to a global nfinimum dependent on the initial guess. If a different step size is used for each parameter, dependence on the initial guess is reduced. Once optimum step sizes have been estimated, new estimates for the coefficients of the parameters in (23)and the leaky parameter in 126) may be computed from
~alue of 126) for the kth iteration. When iteration k is equal to an integer multiple of the total number of parameters being identified, the steepest descent method is applied. This pacer method insures that the parameter identification algorithm does not diverge from the global minimum (Bazarra and Shetty ~1. After the gradient vectors (29) are determined, the line search technique is again applied to locate the optimum step size for each parameter group as before (Bazarra and Shett,v ~ and Gill, et al.~e). Our results confirm the conclusions given by' Neuman :~ and Carrera and Neuman 5 that the conjugate gradient method of Fletcher and Reeves improves the convergence rate of the parameter identification algorithm. Searching for a global minimum continues until the following criteria are satisfied: 1. The maximum number of iterations specified by the user is reached. 2. jk(p) <~ 10 s 3. iAJkip) J k ~(p):! ~ I0 s for 4 consecutive iterations. 5. M O D E L V E R I F I C A T I O N U S I N G S Y N T H E T I C DATA Synthetic data with known parameter distributions with and without random noise were used to test the algorithm. Synthetic data tests the ability of a parameter identification algorithm to identify a known distributed parameter system (Carrera and NeumanS). These synthetic tests may be grouped into the following categories: (1) (2) {3) (4) (5)
Identification of analytical parameter distributions; Stability of the parameter identification algorithm; Uniqueness of the parameter estimates; Effects of scaling on the identification algorithm; Comparison between the steepest descent and Fletcher-Reeves gradient methods; (6) Effects of noise on the identification algorithm.
5.1 hlent(fication of analytical parameter distributions "Measured' piezometric heads h"(r, t) for the two cases
,'= Ii! +::!)I!i where ~.~ is the optimum step size and d~.~is the directional vector for the kth iteration for parameter (-)[Bazarra and Shetty~).
4.3 Fletcher-Reeves gradient method The Fletcher-Reeves rnethod is based on conjugate gradients for a quadratic equation (Bazarra and Shetty'). The steepest descent method generates the first iterate value. For subsequent iterates, a positive multiple of the search direction G~.j from the previous k - 1 iteration is applied to the directional vector for the kth iteration according to
analyzed were computed numerically by the transformed BVP (17). Dimensions and analytical parameter distributions used for the two cases are listed in Table 1. The number of coef'ficients used to estimate the transmissivity for Case I listed in Table 1 by (23a) was N, = 2 and storativity was specified. The number of coefficients used to estimate transmissivity for Case2, Table 1 by (23a) was .,\'~ = 4 and storativity by (23b) was N, = 2. The following Root Mean Square (-RMS) errors ac~ were used to quantify' errors for hlr, t), T{r), and S(r): .
. 1 . x . 1 ~t
~ = ~
m ~] -- hi..)-
[4o
n=l
,
[(h~
L ~
(30a)
i=l
,L [ITI'-TT")-']
a~-= 7 2-
"F-
t30b)
i=I
where
, 'I k
G~.) - -
I2
, -~[ ~~l.~ - - 1 1 ,, 72
I ,L [ ( s [ - s D - ' ]
j~.~ .. o~: . . -.
I.I
(29a,b)
where I(')iJ is a Euclidean norm o f ( ' ) and J~.~ is the
(30c)
i=I
where F/° is the average value of the initial head data H°{r) in (6a); T and S are the average values of Tim and S'~", respectively: I = 25 is the total number of simulated
Adv. Water Resources, 1991. I'~d. 14, No. 5
245
P~lrameter Me, t(1ication in mdial f l o w R. T. Hudspeth et al. 30
25,
25
20
2O
]
I
2" 5
0
~0 F.--
5¸
J
" o-~
0
(a)
0
0
2 4 6 R~zdiaI Oistance
8
1o
0
(a)
20
0.5
10 Rodial
I5 2_.0 Distance
2.5
20
15 >,
~5
!
to ~n
c t3
21,22,25,24,25
IO
E c ID
I0 Theoretical
Theoretical Values
5
(b)
l
o
I
z
i
I
,~
I
i
1
6
i
8
5 ~
7
0
0.5
Vol
,
28,29,50
i
Io
Rodic]l Distance
Fig. 2 Transmissivity estimates for Case I, Table 1 with storativitv speciJied. (a) k = 0-4, (b) k = 21-25 I
measured heads; N = 21 is the total number of discrete time values. Results are illustrated in Figs 2 through 4. The solid triangles on these figures represent the location of the "measured' prior data for T(r) only for Fig. 2; and for both T(r) and S(r) for Figs 3 and 4. In these cases and all others examined, the parameters were well-estimated by the identification algorithm. In most cases the maximum number of iterations criterion # 1 . in § 4.3 was reached before either one of the other two criteria in § 4.3. These results did not depend on whether the simple basis functions (23) or more complex representations were used to estimate the parameters. 5.2 Stahilitv o/ the parameter identification algorithm Stability is achieved, in part, through the adjoint BVP (25) by eliminating the need to differentiate the heads with respect to the parameters. The line search procedure
246
Adv. Water Resources, 1991. Vol. 14, No. 5
(b)
t
I
I
1.0 1.5 2_.0 Redial Distance
2_.5
Fig. 3 Tra,smissivity estimates Jar Case 2. Table 1 • (a) k = O~l. (b) k = 21-25 and variable step size were utilized until one of the convergence criteria in § 4.3 was satisfied. Since the analytical values for the distributed parameters listed in Table 1 were given, convergence was usually obtained during the first four iterations. The program would normally terminate at this point having found the global minimum. To test for stability, however, the algorithm was forced to continue for a total of 50 iterations with the line search procedure omitted and with a step size ~k ~(.~ = 0.1 in (27). Figure 5 compares the initial iterations and final five iterations for transmissivity only for Case 1, Table 1 with storativity specified. The distributed parameter estimates T~(r) in (23a) oscillate about the true analytical value for
Pmametcr lljdllg!/i('d[[¢~ll in radia/jlo~ l R. TI Hu~gpeth et al. 20
J
40
15
>. 30 F. 2
I
>
I, 2,3,4
20
(.,q
c
eorelica[
Values
I0
al Values o.,
0
(a)
I
t
0.5
0
1.0
Radial
z
i
1.5
2.0
2_.5
Distance
50
4o
(a)
0 0
1
i
2
Theoretical Values 30 >.
,
•
t
~
i
4 6 8 Rodiol Distance
I0
£0
27, 28,29,3C
L; *6 20
2 03
15 I0 E c O
l
o (b)
0.5
l
,.o
Radial
i
t
1.5
2.o
a
e.5
~
I0
I..-
"-- 4 6 , 4 8 , 50
,Theoretical Values
Ois~'ance
F(g. 4 Storativitv estimates f o r Case 2, Tahle /." (a) k =0-4. /h) k = 2 / - 2 5
Case 1. The magnitude of these oscillations does not increase as the number of iterations increases proving that the identification algorithm is stable. The magnitude ofthese oscillations is a function of the step size ~.~ chosen in (27).
5.3 Uniquenes.v of parameter estimates Introducing prior data T"(q 0 and S"(~t) contributes to unique parameter estimates. The effects of prior data for only T"(qt) with S(r) specified were examined by the following tests listed in Table 2: Test A and B: Case 2, Table I with S(r) specified and Dirichlet boundary conditions (2 = ~ ) and prior data for only T"bh) simulated at location q~ = 0.3, 762, and 1,524 meters. Test C: Case 2, Table 1 with S(r) specified and Dirichlet boundary conditions {,:. = :c) but only one value of prior data T"Oh) was simulated at qt = 0.3 m at the well boundary.
01
0 (b)
i
i
2
L
i
4 Radial
n
~
L
6 Distance
l
8
i
Io
Fig. 5 Stability oJ transmissivitv estimatesJbr Case 1. Table 1 with storativitv specifi:ed: (a) k = 0-4, (h) k = 46-50 Tests A and B compare the sensitivity of the identification algorithm to initial guesses for parameter estimates. Test C illustrates the need to have prior data T"(th) specified away from the well boundary. Storativity S(r) was specified for all three tests listed in Table 2. The identification algorithm closely approximated the analytical parameter distribution for Case 2, Table 1 for T(r) in both Test A and B independent of the initial guesses. However, the identification algorithm failed to approximate the analytical distribution listed for Case 2,
Adv. Water Reso,trces. 1991, Vo/. 14, No. 5
247
Parameter ideptti/ication in radial flow. R. T HudYpeth et al. TABLE 1. Summary of Test Cases Analyzed CASE i
I 1
i
k
/~S
/2/°
(m)
(m)
(m)
(m~
(m~/d)
0.I52
15,24
152
270
152g
:"
/-
Ss
(:o'2
(m~Id)
~.0 4o+1oo(
2
I
0.152
1524
609
270
(m~/d)
2.0
1.0
P/R
)
2.0
2196
8.0
40+ exp{4.6ff/R)}
I+9( ://~ )
TABLE 2. Summaryof Uniqueness Tests for T(r) only with S(r) Specified (CASE 2, TABLE I) prior da~ ]ocadon:
TEST
l ]
Itc~tion k = 0
Iteration k = 30
(m)
i
i
~,
J(P)
,~,
I
.3
1
I0
I
5.14
2
762
2
0
2
-.97
A
J(p)
361 3
B
C
1524
3
0
3
.05
4
0
4
1.75
1
0.3
l
1.0
2
762
2
0
3
1524
3
0.3
I
RMS ERROR
~l
1
5.06
2
-.57
0
3
.17
4
0
4
1.58
5622
I
1.0
1
5.03
2
0
2
-1.99
3
0
3
-.27
4
0
4
-.17
5476
~k
qr
(~)
(~)
0.61
.2
l.g
.63
,0~
1.7
.40
.I3
27.9
TABLE 3. Effects of Scaling on Objective Function Terms (Case I, Table 1) TEST
(m)
(mild)
A
1.0
1.0
[
1.0
(m)
1.0
(d)
1.0
B
270.
40.
I
I0 "i
1524
2.0
]tcr'zdon k
-:(p)
J,(p)
I
I
I
I
J~)
B
A
B
A
B
A
B
o
5XlO J
3.09
5XI~
2XIO "~
5XIO 1
3
23
-0
A
1
4X10~
0.6
4XI~
6X10 "~
265
.62
.16
2
472
-0
4-45
-0
27
-0
2XIO 2
3
.27
.25
.01
1XIO 2
4
3XIO ~
2XlO"
lXl04
2XI0'
5
18XI0 "s
2XI0 ~
6XI0 7
2XI0 "~
6
16XIO s
2XIOa
6XI0 "7
--0
7
25XI0 -j
3XIO ~
6Xlff ~
2XI0 ~
6XI0 "~
8
18XIO"#
9
18XIOq
10
10XI0 "I
I"
6XlO a
2XI0 ~ LXI0 ~
Table 1 for T(r) in Test C when only one value for T"(r/t) was provided near the well b o u n d a r y (qt = 0.3 m). The parameter distribution for T(r) for Case 2, Table 1, estimated by Test C, Table 2 is illustrated in Fig. 6. N o t e that the final value of the objective function j3O(p) in Table 2 was nearly the same for all three tests; but the RMS errors for a T were not. Test C clearly demonstrates
248
:,(p)
Adv. Water Resources, 1991, Vol. 14, No. 5
-
6XI0 a
the inability of the identification algorithm to find unique parameter estimates without prior data T'(qt) being specified away from the boundaries when Dirichlet b o u n d a r y conditions (). = z ) are used. Again, storativitv S(r) was specified for all three cases listed in Table 2.
5.4 Effects of scaling on the identification algorithm
Parameter ident~ication in radial flo~ 20
io,ooo
-
R. T. Hu~£peth et aI.
-
I,O00
15 =x
g
>
E
er-Reeves
t~ c
I0
cal c o
Ioo
Z Theorelicol
Values
on
>
Io
u
',7 0
I
26,27, 2 8 , 2 9 , 3 ~ ~
0
o
(a)
I
i
0.5
I.o
Radial
0.I
i
2.0
1.5
Distance
i__
o
Io
20
i
30
L.__..
40
50
r terotions
Fig. 7 Comparison between steel)est Fletcher- Ree yes ,gradient methods
descent
and
20 ten iterations listed in Table 3. With the scaling used in Test B, Table 3, the error in head values Jh(P) was always insignificant compared to the prior data term Ja(p). For Test B. the indentification algorithm converged in just two iterations. However, the head errors Jh(P) were never greater than 10 -s because of scaling and contributed essentially nothing to the objective function (18). Identification of the transmissivity parameter for Test B was entirely the result of the prior data term Ja(P) because of scaling.
15
E I0 • Theoretical Value
I--
5--
5.5 Comparison between the steepest descent and the Fletcher- Reeves gradient methods
4../-
o-N i
0
(b)
0.5 Radial
I
1.0
1
1.5
,
2.0
Distance
Fig. 6 Uniqueness Jar Case2, Table 1 b3' Test C, Table 2. (a) k = 0--4, (h) k = 26-30
The effect of scaling on the rate of convergence may be significant. Two related scaling problems were evaluated: (1) the effect of improper scaling; and (2) the effect of scaling the basis functions. Data for Case 1, Table 1 with S(r) specified and Dirichlet boundary conditions (/. = m) were simulated for the two tests listed in Table 3. Test A has the effect of not scaling the algorithm; while Test B scales h~.,, S~, R~ and t, by their largest values and T~ by the reference value at the well. These two tests represent two extremes which could be applied. For both tests, prior data values for only T'(q~) were again specified at locations q~ = 0.3, 762, and 1,524 meters. The effects of not scaling terms in the objective function (18) are illustrated by Test A, Table 3 where the error in head values Jh(P) is two orders of magnitude larger than the prior data term Ja(P) for the first iteration k = 0. This difference in magnitude remained constant for all of the
Using Case 1, Table 1 with S(r) specified and Dirichlet boundary conditions (2 = :c) as input, the identification algorithm was run twice using the steepest descent method first and the Fletcher-Reeves gradient method second. The results of this comparison are illustrated in Fig. 7. The steepest descent method converged to a global optimum fastest.
5.6 Effects of noise on the identiJicatio,z algorithm Effects in measured heads hi.",, were simulated by introducing Gaussian white noise into the synthetic head data only. The identification algorithm searched only for the transmissivity distribution T(r) for Case 1, Table 1 with S(r) specified. Prior data values for transmissivity T,"~ were simulated without noise at locations rh = 0.3, 762, and 1524 meters. To evaluate whether the noise was truly random, the RMS error and maximum error in the 'measured' heads was evaluated using the theoretical values for Case 1, Table 1. The five tests listed in Table 4 simulate increasing amounts of noise. If the head noise is Gaussian, then the resulting RMS error a~ listed in Table 4 should be approximately 50 % of the maximum noise. For examples listed in Table 4, these conditions were satisfied. Table 4 indicates for Case 1, Table 1 with :~, = 5.0 and :~_, = 2.5 with head noise that the algorithm could identify the transmissivity as large as + 0.152 meters (Tests A-C).
Adv. Water Resources, 1991, Vol. 14, No. 5
249
P a r a m u w r ,dumiHuation in r¢Mial f l o w . R. T. HudYputh et al. TABLE 4. Effects of Noisy Head Data on T(r) Specified (Case I, Table i) M rasu rea:l Hem.ds
Itermtaon k = 0
~loise
Teat
0"n)
(%)
d (%)
t
Me, x E r r o r
i
~,.I
[
i%)
crr
o,
t
Ii
A
5.0xl0 ~
±0.003
6.0x10 a
-2.@,lxl0
6.88x10 a
.5.53x10 ~
0.016
a..99 2'.',9
3.09x.tO "~
15
B
5,Ox10 3
±0.030
5.3x10 ~
-2.6ax10 2
6.74x.I 0 ~
4.8ax.I 0"~
0,O65
4.99 2.50
0.289 t
20
C
2.Yx10 "i
+0.152
0.02~,
41.135
0.0284
4). 139
0.660
500 2 45
6.995
25
D
0.05
±0.305
0.056
0.269
0.136
1.02'7
1.058
5.12 2.50
30.a2
25
E
0.28
± 1520
0.312
1.50'0
0.,196
4.06~
2.730
aL .65
1055.8
25
,I
2.60
I
F o r Test E with head noise ___1.52 meters, the a v e r a g e p e r c e n t e r r o r in the t r a n s m i s s i v i t y e s t i m a t e s after 25 i t e r a t i o n s was only %. = 3 % .
REFERENCES I
5.7 Com'ht,vion,~
2
1. T h e p a r a m e t e r i d e n t i f i c a t i o n a l g o r i t h m successfully a p p r o x i m a t e d the t h e o r e t i c a l p a r a m e t e r d i s t r i b u t i o n s .
3
2. T h e p a r a m e t e r i d e n t i f i c a t i o n a l g o r i t h m was stable with small, b o u n d e d o s c i l l a t i o n s a b o u t the true t h e o r e t i cal solution. Stability was o b t a i n e d by the use of the a d j o i n t BVP. 3. T h e i m p o r t a n c e of p r i o r d a t a was illustrated in Fig. 6 w h e r e the correct p a r a m e t e r d i s t r i b u t i o n was not identified e v e n t h o u g h the o b j e c t i v e f u n c t i o n was m i n i m i z e d . F o r this case, prior d a t a was r e q u i r e d in o r d e r to o b t a i n a uniqt, e solution.
4 5 6 7
4. T h e effects of scaling on the o b j e c t i v e f u n c t i o n ilhtstrated that i m p r o p e r scaling results in s o m e t e r m s h a v i n g no c o n t r i b u t i o n s to the o b j e c t i v e function. P r o p e r scaling insures that all t e r m s in the o b j e c t i v e f u n c t i o n are comparable.
9
5. B o t h the steepest d e s c e n t and F l e t c h e r - R e e v e s gradient m e t h o d s c o n v e r g e d to a g l o b a l m i n i m u m .
10
6. T h e results for noisy head d a t a illustrated that the p a r a m e t e r identification a l g o r i t h m is stable and will c o n v e r g e to the t h e o r e t i c a l p a r a m e t e r d i s t r i b u t i o n e v e n with r a n d o m errors in the m e a s u r e d heads. T h e m a g n i t u d e of the o b j e c t i v e f u n c t i o n was n o t a g o o d m e a s u r e of the g o o d n e s s of fit for the p a r a m e t e r s . T h i s w o u l d suggest that several trials are r e q u i r e d for a p p l i c a t i o n s to field d a t a w h e r e the d i m e n s i o n a l i t y of the basis f u n c t i o n s and the radial scale m u s t be a d j u s t e d in o r d e r to d e t e r m i n e an o p t i m u m result.
ACKNOWLEDGEMENTS W e gratefully a c k n o w l e d g e the financial s u p p o r t p r o vided by the D e p a r t m e n t of E n e r g y ( D O E ) u n d e r the Basalt W a s t e I s o l a t i o n P r o j e c t ( B W I P ) to the N o r t h w e s t C o l l e g e and U n i v e r s i t y A s s o c i a t i o n for Science ( N O R C U S ) u n d e r c o n t r a c t D E - A M 0 6 - 7 6 - R L D - 2 2 2 J ; a n d the Office of N a v a l R e s e a r c h ( O N R ) u n d e r the U n i v e r s i t y R e s e a r c h Initiative ( U R I ) C o n t r a c t N 0 0 0 1 4 - 8 4 - C - 0 3 0 3 . W e h a v e benefitted f r o m n u m e r o u s s t i m u l a t i n g discussions with Drs. S t e v e B a k e r a n d Allen Liu.
250
M a x Err,or (m)
Adv. Water Resources. 1991, Vol. I4. No. 5
8
I1 12 13 14 15 16 17 18 19 20
Bazarra, N1. {5. and Shctt}, C. M., :\mllinuur Prowammin~l, Theory. and -II¢torirhm~. Wile}}} and Sons Inc., New York, 251 33(I. 1979 Carnahan, B., Luther, H. A. and Wilkes, J. O.. Applied :V,mwrical 3,1ethod,~, John Wile}}}' and Sons, New York, 113 115 and 429 461, 1965 Carrcra. J. and Ncuman, S. P, Estimation of aquifer paran'~etcrs under transient and steady state conditions: I. maximum likelihood method incorporating prior infc,rmation, II'RR, 22. 199 210, 1986 Carrera, J. and Ncuman, S. P. [:.stimation of aquifer parameters under transient and steady state conditions: 2. uniqueness. stability, and solution algorithms. WRR, 22, 211-227. 1986 Carrera. J.and Neuman. S. P. Estimation of aquifer parameters under transient and stead,, state conditions: 3, application of synthetic and lield data, WRR, 22, 228 242. 1986 Chavcnt, O. History matching by use of optimal theor,',, S. 0! Pet. Knelt.. 74 86. 1971 Cochran. J. C., Wiser. W. C. and Rice, B. J. ,-tdl'anced Enqineering Muthemuti{s, Brooks,Cole Publishing, Monterey. CA, 175-184, 1987 Cooley, R. l_. Incorporation of prior information on parameters into nonlinear regression groundx~ater tlov. models: 1. theor,,,. WRR, 18, 965 976, 1982 Ehlig, C. and Halepaska, J. C., A numerical study of confinedunconlined aquifers including etlkzcts of delayed yield and leakage. WRR. 12. 1175 1183, 1976 Friedman. B. Principles and Techniques Ol Applied Mathematics, Wile,,' and Sons, New York, 42 147, 1956 Gerlach, J. and Guenther. R. Remarks on part, meter identification, Convergence of the algorithm, Numerical .%lathemutics, 51, 3 9, 1987 Gill, P. E., Murray, W. and Wright. M. H. Practical Opfimizution, Academic Press. 1986 Guenther, R. B., Hudspeth, R. T.. McDougal, W., Gerlach, J. Remarks on parameter identification, I. The algorithm. Numer. 3,lath. 47. 355-361, 1985 Guenther, R. B. and Lee, J. W. Partial D!fferemial Equations of .~,lathematical Physics and Intewal Equations, Prentice-Hall, Englc;vood Cliffs. New Jersey. 16, 19 and 362, 1988 Hadamard, J., Le Prohl&ne de Cauehy et les Equations aux Dt;rivdes Partielh's Lineaires ffyperholiques, Hermann, Paris. 1932 Hildebrand. F. B. Finite-Dffl~,rence Equations and Sim,dutions. Prentice-Hall, Inc., p. 221. 1968 John, F. Partial differential equations, Ch, V, Mathematics Applied to Physics. Springer-Verlag. Ne,a York, NY, 273 277, 1970 Kravaris, C. and Seinfeld, J. Identification of parameters in distributed parameter systems by' regularization. SIAM J. Control and Optimization, 23,217 291, 1985 McWhorter, D. and Sunada. D. K. Ground~tater ttydrology and Hydraulics. Water Resources Pub.. Littleton, Co. 177-194, 1977 Neuman, S. P. Calibration ofdistributed parameter groundwater tto~v models viewed as multiple objective decision process under
Parameter identification in radial J t o w
21
22
23
24
25 26 27
uncertaint), WRR, 9, 1006-1021, 1973 Neuman, S P. A statistical approach to thein,,erse problem of aquifer h?drolog?: 3. improved solution method and added perspecti,,e, WRR. 16. 331-346, 1980 Neuman. S. P. and Yakowitz, S. A statistical approach to the in;erse problem of aquifer h.,,drology: 1. theor,,, ~IRR, 15, 845 860. 1979 Neuman. S. P., Fogg, G. E. and Jacobson. E. A. A statistical approach to the inverse problem of aquifer hydrology: 2. case stud,.. WRR, 16, 33-58, 1980 Role). K. L. An indirect parameter identification algorithm in radial coordinates for a porous medium, MS 77wsis. Dept. of Cil'il Enqmeerin~l, Oregon State University, Corvallis, O R. 1992 Ruston. K. R. and Redshaw, S. C. Seepa~le aml Grolmdwater Flm~, Wile~, Nev, York, NY. 32 35. 1979 Stakgold, l, Green's Functions and Boundary V~due Problems, John Wiley & Sons. New York, NY, 43, 1979 Sun. N. and Yeh, W. W. Identification of parameter structure m groundwater inverse problem, WRR, 21,869-883. 1985
28
29
30 31 32 33 34
R. 7-. Hudspeth et al.
Theis, C. V., The relation betv,een the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, Tr~ms. .4mer Geophys. Union. 2, 519-524, 1935 Tikhonov, A. N. Solution of incorrectly formulated problems and the regularization method, Sot. Math. Ookl., 4, 1035-1038, 1963 Tikhonov, A. N. Regularization of incorrectly posed problems, Sot'. Math. Dokl., 4, 1624-1627, I963 Vemuri, V. and Karplus, W. J. Digital Computer Treatment of Partial D(~'erential Equations, Prentice Hall. Englewood Cliffs, NJ, 281-282, 981 Willis, R. and Yeh, W-G. Groundwater Systems Plannin9 und Munu~lement, Prentice Hall Publishing, 1987 Yakowitz, S. and Duckstein, L. Instability in aquifer identification: theory and case studies, WRR. 16, 1045-1064, 1980 Yeh, W. W-G. Review of parameter identification procedures in groundwater hydrolog,', : the inverse problem, WRR, 22, 95-108, 1986
Adv. Water Resources, 1991, Vol. 14, No. 5
251
Department oj" Civil Engineering. Oregon State University. Corvallis. Oregon 97331. USA The classical inverse parameter identification problem is solved using an indirect method. This method is based on the minimization of an objective function or error criterion consisting of three parts: 1) least-squares error of head residuals; 2) prior information of flow parameters: and 31 regularization. An adjoint equation is incorporated into the method to eliminate the need to differentiate the heads with respect to the parameters being identified, increasing the stability of the algorithm. This method is applied to a radial coordinate system used for describing an anisotropic confined aquifer. A variable implicit finite difference method is used for predicting hydraulic heads given a continuous set of flow parameters. The method is tested using synthetic data and is found to give estimates which are both stable and unique.
1. I N T R O D U C T I O N
using only values of head according to
During the past 25 years, a number of algorithms have been developed to solve the classical inverse problem for groundwater flow. Yeh 3"*and Willis and Yeh 32 reviewed these algorithms in detail. Neuman e° has proposed classifying parameter identification algorithms as either direct or indirect. The direct algorithm is a non-iterative technique where the parameters are unknowns and the head values are knowns. Since piezometric head measurements are usually available only at a few specified locations in the field, some means of interpolation must be used to estimate heads at other spatial locations resulting in interpolation errors. The results obtained by the direct algorithm also seem to be correlated to the grid size used {Willis and Yeh32). The indirect algorithm uses an iterative technique to minimize an objective function. The least squares error between the measured and estimated heads provides the basis for many objective functions (Neumana'; Cooley~; Sun and Yeh-'r; Carrera and Neuman3"*'s). Additional penalty terms mav be added to the least squares error resulting in a multiple objective function. The least squares criterion is capable of filtering out some of the measurement errors inherent in field data. However, this term alone is not enough to guarantee a unique set of parameters (Carrera and Neuman3). A unique set may be computed by the addition of a prior data penalty term to the objective function. The inverse parameter identification problem may be conceptualized as follows: Given the fvnctional relationship
h(r, t ) = For [p(r)]
(1)
where h(r. t) is the hydraulic head; r and t are space and time, respectively; and plr) are spatially varying parameters that control the groundwater flow: identify p(r) * Department of Mathematics, Oregon State University, Corvallis, Oregon 97331, USA C 1991 Elsevier Science Publishers Ltd.
240
Adv. Water Resources, 1991, Vol. 14, No. 5
p(r) = Inv [h(r, t)]
(2)
Carrera and Neuman'* identify (1) as a "forward relationship' and (2) as an "inverse relationship.' This inverse problem is well-posed in the sense of Hadamard 15 if: (1) For every value of head data h(r, t), there is one and only one solution for the parameters p(r), and (2) The parameters p(r) depend continuously on the head data h(r, t) (John 17, Guenther and Leer'*). Condition (1) insures that the parameters are identifiable and unique. Identifiability is not achieved in the forward relationship if different sets of flow parameters result in the same head distribution. Uniqueness is not obtained in the inverse relationship when different head distributions produce the same set of flow parameters. This distinction is important since an indirect parameter identification algorithm may produce an identifiable solution that is non-unique (Carrera and Neuman'*). Condition (2) insures that the values of the parameters are stable. The parameter values are unstable if small errors in the head data result in large errors in the parameter values. Two types of errors are known to cause instability: 1) modeling and interpolation errors, and 2) measurement errors (Neuman and Yakowitzee). Analyses of the sensitivity coefficients indicate that errors in the head data contribute the most to instability. The sensitivity coefficients are defined as the change in head with respect to the change in the parameters being identified. Values of sensitivity coefficients that are nearly zero indicate for the inverse relationship that a small change in measured head (i.e. measurement error) may result in a large change in the parameter values (Yakowitz and Duckstein33).
1.1 htentifiabilit.v Sufficient values of measured head data are assumed to be available so that a measurable head gradient exists between any two spatial locations within the flow region during a finite period of time for the duration of the test (Carrera and Neuman'*). It is not always possible to satisfy
Parameter identf/ication in radial flow . Q~(~)
R. T. Hudspeth et al.
.*P~ezomefer
~.-:..:.~-~.: :."..:~.,.q ~ : ~ , T X ~ Z ,
j [:.,.....c..;:.,_..,'.'..'.'.,.
/
".
r=
"
:h =
:t =
" Sir)=
• T(r)=~(3a-e)
E
[\
/
h° {r.lLll
'/
h(r.l~.l) I ~".~
.~," ".7
L(r. ll.l) "\\~'\~x
Flow
Oomo,n
'\ " "
"\
i
The parabolic diffusion equation governing radial flow through porous media shown in Fig. I is well-known (McWhorter and Sunada ~9) and is given in dimensionless form by ~ q h ( r , t)) =
~N'.
F~,. 1.
\ \\
S(r)
~h(r, t)
?t
1 ~ ( (:h(r, t)) Pr ,, rT(r) ~ = 0; ('r
r
\"
r" < r < R . t > ~ O
Dqfinitio~z sketch
this assumption on the boundary of the flow region: and thus the identifiability of the parameters cannot always be guaranteed. Dimensionality may be reduced by estimating the "true" parameter distribution with functions having finite number of coefficients (Carrera and Neuman4). 1.2 Uniqueness
(4)
in which ~'(. ) is a linear operator, S(r) is storativity, T(r) is the transmissivity, r" is the well radius, and R is the radial distance to the boundary of the flow domain. The initial and boundary conditions may be scaled by riO(r) /~h(t) h°(r) = Ft~ " H a ( t ) - f'f (~'(0. Qh(t)" Q'(t) = ~ , Qb(t) = ~ ;.
)~ T~
(5a-e)
The dimensionless initial and boundary conditions are
Unique parameter values result from the indirect method when prior data are included in the error criterion (Neuman2°). Neuman et a[. 23, Guenther, eta[. ~3, Carrera and Neuman 4's and Gerlach and Guenther ~ demonstrate the use of prior data. Prior data are parameter values that have been "measured" or estimated at specific locations within the flow domain. Prior data may be obtained from eflher field or lab tests or by comparison with similar aquifers. It must be assumed that all measured parameter values contain some noise, the magnitude of which may or may not be known. 1.3 Stability Stability may be obtained by several methods. The adjoint equation (Chavent 6) may be used to eliminate all sensitivity coefficients from the gradient equations used to minimize the objective function. The relationship between instability and the sensitivity coefficients is demonstrated by Yakowitz and Duckstein 33. Some algorithms that use the adjoint equation do not eliminate the sensitivity coefficients from the gradient equations used to minimize the objective function (Sun and Yeh2:). An additional penalty term composed of the head residuals at the final value of time may be added to the error criterion. The addition of this penalty term improves stability by taking advantage of more reliable values of measured final heads due to repeated measurements and the reduction in the rate at which the head values are falling. Finally, Kravaris and Seinfeld ~s and Gerlach and Guenther ~ add a regularization term to the objective function to improve stability.
h(r, 0) = H°(r), lira
rT(r)
~h(r, t) ?r
rT(r) -?h(r, t) ~r
-
r"'< r < R
(6a)
-Q'(O
(6b)
'
t > 0
2~
,;. 1 27z h(r, t) = ~ [--,;.Hb(t) -- Qh(t)]; r = R, t > O
(6c)
where tfh(t) is the prescribed head at the radial boundary; Qh(t) is the prescribed flux across the radial boundary; and ). is a leakage coefficient at the radial boundary. Note that the leaky, head-driven boundary condition given by (6c) may be used to prescribe Dirichlet boundary conditions when ). = :,-. and Neumann boundary conditions when ). = 0. Numerical solutions may be improved if the arbitrary time-dependent pumping function Q'(t) is defined by N
Q"(t) = y" [Q~ + Qii(t)].U,it. r, ,, r,,)
(7)
n=I
where N is the total number of finite time intervals during the total test duration; Q~, is a constant pumping rate during time interval n; Q~(t) is a variable pumping rate during time interval n: and the Heaviside Step function U,(t, ~,,_ ~, L,) is a pathological numerical switch defined as (Cochran et a/.7), U,(t, v,_ 1, ~,) = U(t - z,_ 1 ) - U(t - z,)
(8)
that turns the well pumping function Q'(t) on at t = v,_ and off at t = % during the time interval n. 2. la Hybrid analytical-numerical solution
2. B O U N D A R Y
VALUE
PROBLEM
(BVP)
Dimensional variables denoted by (?) are scaled by dimensional scaling parameters denoted by (')s to obtain dimensionless variables according to
The strategy used in the hybrid analytical-numerical approach is based on linear superposition (Stakgold26). The solution to the dimensionless BVP (4 - 6) is assumed to be a linear combination of the following linearly independent functions:
Adv. Water Resources, 1991. Vol. 14, No. 5
241
Parclme[cr idcnHticatiim tn rcldial jlon. R. T. HuJ~pcth et al. hit. ti = h ] [ I C ] :
[BC~]: [BC:]: [I-]', _
= h,,:[IC]:
0: 0: 0', + hi',0: [BC,]:
0: 0',
where [IC] are initial conditions: [-BCi] are boundary data at the ith boundary': and [ j ] is the forcing function. The linear superposition principle suggests looking for a linear combination of linearly independent functions which satisfy only specific elements of the inhomogeneous system while contributing nothing to the other inhomogeneous elements (i.e. they' satisfy homogeneous conditions elsewhere). In principle, a linear superposition approach would seek the linear combination (9b): however, in order to take adwmtage of the well-known Theis integral (Theis es, Guenther and Lee~'~), the ideal linear combination (9b) will be replaced by' (I0)
where h°(r, r) satisfies the initial conditions: and ,4r, t) and r(r, t) are analytical and nurnerical approximations to the b o u n d a r y data. Substituting (10)into (4) gives 5/'(h) = 5/'(h °) + 5/'(u) + 5/'(r) = 0
(l I)
The initial data are satisfied exactly by the solution to 5/'(h°): while the b o u n d a r y data are satisfied by the solutions to 5/'(u) and L/'(v). The following dimensionless BVP is defined for t/~(r, t): L/'(IC) = S(r) Fh°(r, t)
I C IrT(r) (h°(r, = t) 1
(:t
=0:
r" < r <
r 7r
(12a)
R,t>0
with inhomogeneous initial data h°(r. 0)= H°b'):
r'"
R
(12b)
and homogeneous boundary data ~:h°(r,
t
0:
c~r
r= R:t >0
(12c)
The following dimensionless BVP with constant coefficients is defined for u(r, t):
~,O(u)=So&4r.
t)
(':[
T° F (Cu(,', t)) -
=
I"
r
.
(T
( F
=0:
,
r'" < r < R, t > 0
(13a)
with homogeneous initial data u(r, 0 ) = 0,
r " ~ < r ~< R:
(13b)
with inhomogeneot, s b o u n d a r y data at the well lim r
r' ~ 0
rT o ?u(r, t) dr
-
-Q"(t) 2=
t >
O:
(13c)
and with a homogeneous far field condition u(r,t)=0.,
r> R,t >0
(13d)
The parameters S o and T o are values of these parameters at the well r = r" --+ 0. The dimensionless solution to BVP (13) is the well-known Theis integral given by
242
Adv. Water Resources. /991. Vol. 14. No. 5
t't
o O- " ( t ' )
S(r)
Ft'(r, t) ~t
l ((rT( v r- ~r
(14)
(r - t)
&'(r, t)) r) ~ : T f - - ,
. =
l(r,
t)
r" < r < R , t > 0
(15a)
t)] -~:u(r" - - - t) + I .F ("rT(r)(U(r. . . .
(15b)
with forcing function J(r,t)=-S(r)
~t
r ~r
~'r /"
with h o m o g e n e o u s initial condition dr, 0 ) = 0
r" ~ r ~ < R
(15c)
with inhomogeneous b o u n d a r y data at the well lira r=r"
rT(r)
FRO', t)
-
~h°(r, t) rT(r) - -
t > 0
(15d)
~0
and with inhomogeneous far field b o u n d a r y data ~:v{r,t)
rW(r)
Fr
}.
1
2= v(r,t) = ~ - rT(r)
[
}.Hh(t) -- Qh(t)]
~:u(r, t) " r = R. t > 0 ~r
(15e)
The solution for h(r, t) is now defined by the following linear system of dimensionless BVPs: h(r, t) h°(r, O)
{[IC]: [BC,]: EGG] [0]I
~r
2
The dimensionless BVP for r(r, t) becomes
(9a,b)
h(r, t) = ttr~(r, t) + u(r, t) + dr, t)
i
ulr, t) - 4]4-7%
+ hj.',O: 0: 0: [./]i
It_,',(): 0: [BC:]:0;
-
Sot
+
u(r, t)
',[o]: [BC, ]: [0]: [0]',
+
',E/c]: [o]: [o]. [o]', r(r, t)
',[o]: [BC,]; [~C_,]; [.I-]I
The Theis integral was c o m p u t e d using both infinite series and G a u s s - L a g u e r r e quadratt, re (Roley2a). The forcing function f(r, t) in { 15b) must be calculated from u(r, t). The value of f(r, t) represents the residual to (ll) when the distributed parameters T(r) and S(r) are not constants. If T{r) and S(r) are constant, then f(r, t) =- O. The hybrid analytical-numerical method offers several advantages over the totally numerical method. (1) Lower sensitivity to At. During the initial interval of the p u m p test, the analytical solution by the Theis integral (14l is independent of At and accounts for most of the drawdown. Therefore, regardless of the p a r a m e t e r distribution, the initial solution to (I 1) is nearly analytical with r(r, t) --- 0. This is a particular advantage because d r a w d o w n is most rapid at the initiation of pumping, and numerical errors tend to be greatest during this interval. (2) Smaller numerical errors. The numerical solution r(r, t) is only' required to correct the Theis integral (141 when the parameters are non-constant or distributed. Sensitivity' to numerical methods is thereby reduced. (3) Separation of variables. The linear superposition given by {16) provides a convenient way to separate the effects on the Theis integral of an aquifer having leaky b o u n d a r y conditions, distributed parameters, etc.
Parameter identification in radial flow. R. 7". Hudspeth et al. A disadvantage of the hybrid analytical-numerical method is an increase in the computational time compared with a fully numerical method. For example, a solution using 75 spatial nodes and 50 time steps on a personal computer required 140% more computer time.
2. lb Numerical method The radial derivative ~h(r, t)/& in (6b) approaches infinity at the well as a result of a singularity at r = 0. To approximate this derivative accurately near the well boundary in discrete form, a very small Ar is necessary. A coordinate transformation that was introduced by Ehlig and Halepaska 9 and used by Ruston and Redshaw 25 maps a variable Ar into a constant At) by
for transient data used by Carrera and Neuman ~'. Both of these objective functions include weighted least squares error and prior data terms. The objective function 118) differs from most other objective functions by the addition of the following terms: (1) A weighted least squares error in the measured and predicted heads, Jflp), at the final value of time, r x, which optimizes the improvements in data measurements made at later times. (2) A regularization term, Jg(p), which smooths the parameter estimates if errors are present in the measured data.
The dimensionless BVP (4 - 6) in r is transformed to
In addition to adding these two terms to the objective function, the BVP ( 4 - 6) is scaled in order that each error or parameter estimate in the objective function are of the same order of magnitude. This results in a faster convergence rate for the parameter identification algorithm.
(r"2exp[2p]S(P) g'h(p'
3.1 Least squares error, JnfP)
[)p~ [T(p) ch(p' t)] = O
(17a)
A weighted least squares error between the "measured' or field heads, h"(r, t), and the predicted heads, hP(r, t), is defined by
{17b)
Jh{P) = ~
with dimensionless initial conditions
h(p,O)= H°(p);
0
P
lira T(pJ ,,=~
T(p)
?h{p, t) ?p
~h(p. t) (p
- Q"(t) - - - " 2r~ '
Vl(r, t)[hP(r, t) - h (r,t)]-rdr;
dt L
and dimensionless boundary conditions
r"
1 ~n<~N t>0
(17c)
) h(p, t)= 1 27z ~ [-- )'Hb[t) - Qb{t)]; p= P,t>0
(17d)
where P = In (R/r'") and ?h(p, t)/@ at p = 0 is now finite. The dimensionless transformed BVP (17) may now be solved very etticiently by numerical methods.
2. lc Numerical solution A numerical solution was computed using a weighted implicit finite-difference method (Carnahan, et al.2). The value of the weighting factor 0 may be chosen such that when 0 = 0, a fully explicit method results; when 0 = 0.5, the Crank-Nicholson method results; and when 0 = 1, a fully implicit method results. Numerical results indicated that 0 = 0.6 provided the best estimates. At p = P, a fictitious node is defined outside the domain using the boundary simulation technique described by Hildebrand ~6. The boundary condition at p = P is then used to eliminate the fictitious node. The same technique is also applied at the well, p = 0 (Roley-"~). The dimensionless BVP (17J is reduced to a set of algebraic equations in which the coefficient matrix is tridiagonal. The Thomas algorithm (Vemuri and Karplus 3t) was used to solve this matrix. 3. O B J E C T I V E F U N C T I O N The objective function consists of four terms defined as (Carrera and Neuman "t)
J(pt = Jh(P) + Ja(P) + Jf(P) + JR(P)
(18)
The objective function (18)is similar to equation (24)
(19)
where V~(r, t) is a weight which represents the confidence in the measured head data h"(r, t). The least squares error is the basis for most indirect parameter identification algorithms. However, this term alone is not enough to insure either stability or uniqueness in the parameter estimates.
3.2 Prior data error. Ja(P) In order to achieve uniqueness, a prior data error is defined by 1 ~
Jd(P) = ~ L -
,
V2t[Se(~t) -- sm "(~t)]-
/=I
+S
1
V3t[Tt{tlt) - TmOhl] 2
(20)
- l = I
where S"(~ t). . . . . S"(~.,t) and T"(qt) . . . . . T"(qr) are measured values of storativity and transmissivity, respectively, measured at specific locations within the confined aquifer (i.e., ~l . . . . . ~.~t and ql . . . . . qr); and Se(~l) . . . . . Se(~,,t) and Te(q~) . . . . . T~(qr) are the estimated values of storativity and transmissivity, respectively, obtained from a parameter identification algorithm. The weights V2t and V3~ represent the confidence in the measured prior data. The addition of prior data may not reduce the least squares error observed in the head values. Because of possible errors in the measured head data and the difficulty in obtaining accurate estimates of the flow parameters from field data, the minimum value of the objective function may even increase as a result of adding prior data (Neuman2~). Previous attempts to compute unique parameter values based on experimental data and/or professional judgement have yielded reasonable parameter estimates. The two remaining terms in the objective function (18) are added to provide stability.
Adv. Water Resources. 1991, Vol. 14, No. 5
243
Parameter ident(/ication in radial flow. R T. Hmtspeth et al. 3.3 Final head error. Jr, p,,
The heads measure at the final value of time L~ may be more reliable than earlier measured values because of repeated calibration checks and improved data recording procedures. Accordingly, these final measured heads are weighted more heavily. If the reliability of measurements at a particular location r~ are uncertain, the weighing coefficient l:~.(r~)for these locations may be set to zero.
descent and Fletcher-Reeeves gradient methods (Bazarra and Shetty~). First. the adjoint BVP is derived using variational calculus to eliminate the sensitivity coefficients from the gradient equations. Then. the steepest descent and Fletcher-Reeves discrete gradient methods are used to minimize the objective function. Two potential difficulties ma? be encountered when calculating the gradient equations. First, these gradient equations include sensitivity coefficients which may lead to instabilities. Second, these sensitivity coefficients must be calculated from a BVP requiring the solution of N t + N, + l differential equations. The adjoint BVP may be used to eIiminate both of these potential difficulties.
3.4 Tikhom~v regularization. J~(p)
4.1 .q~Ooittt B I'P (Backward Heat Equation)
Since instabilities in the values of parameter estimates appear to be closely' related to noise in the measured data, a Tikhonov regularization term is added that smooths the parameter estimates. Although there are several ways this could be done, the method of Tikhonor :~''~° used by, Kravaris and Seinfeld Is and by Gerlach and Guenther ~1 is used and is defined by
The adjoint BVP may be derived by first multiplying the variational of the state equation (4) by a scalar function w(r, t), and then integrating by parts over both the spatial and temporal domains to obtain (Friedman ~°)
A weighted least squares error between the measured and the predicted heads at the final value of time is defined bv I
" R
J /(pl = ~_..... Va(r)Slrl[hPlr. r~.) -- h"(r, r~.)]2rdr
= 2
+
121)
S,.,,.) ~',,'(,', - - t) + -1 ( f-~. )
?t
Jrdr
(22)
where Vy(r)and l/~,(r)are weights. This regularization term imposes a penalty on oscillations in parameter estimates by increasing the value of the objective function (Carrera and Neuman3l. The use of the adjoint equation provides an important contribution to stability of the parameter identification algorithm. The adjoint equation can be used to eliminate the need to differentiate the measured heads with respect to the parameters being identified. Note that all of the error terms, Jh(p), Ja(P), and J/(p), reduce to zero when the predicted and measured values are identical. However, for the regularization term (22), J~(p) is always positive-definite. This is not a problem since the objective function is minimized relative only to an initial value. This point is important when establishing criteria for terminating iterations in a parameter identification algorithm.
3.5 Reduction of dimensionality of parameters The distributed parameters, T(r) and S(r), may be approximated by constructing sequences using simple polynomial basis functions according to
j=l
,,,,.flrI" I~
4. M I N I M I Z A T I O N FUNCTION
The objective function is minimized by both the steepest
244
Adv. Water Resources. 1991. Vol. 14. No. 5
~r
(25a)
with inhomogeneous initial data w(r, r x ) = -l,;,(r)[hQr, r s) - h " ' ! r , rs)]:
r'~
and homogeneous boundary data lim
rT"lr)
~;w(r, t)
-0:
0
(25c)
2re w(r. t): 0 < t < r x, r = R
(25d)
~2r
r=r"~O
and
rTqr)
~w(r, t)
(r
)."
where ): is a leaky parameter that must also be identified. If V dr) =- 0 in (25b), the homogeneous initial condition used by Carrera and Neuman "~5 is obtained. The solution to the adjoint BVP (25) reduces the first variation of tile objective function (18)to
dt ill ~hP("" = t)(~,c[,-, t)
J(.)(p) =
(T
~~r
dt -R Hip(r, ( ~ - t)
w(r,
j
(23a,b)
OF THE OBJECTIVE
I rTqr) ,'w(,', t)]
r'"~r~R,O
m=l
where the estimates Tqr) -+ T(r) and S"(r)---, S(r) as N i, N_, ~ m. Parameter identification is now reduced to identifying a finite number of coefficients, :~j and tim, in (23) such that the objective function (18) is minimized.
r ~r
= V~(r, t)[/f'(r, t ) - h"(r, t)]:
N.
Tqr)= ~. :~j"r Ij 1~., S~'(r)= 2
(24)
where ~'*(.) is the adjoint operator. The inhomogeneous adjoint BVP for w(r, t) is
w
/~Tqr)~2q
N ~
T~',O')rdr
[)S[.)(rt
rdr
r~
.If
+E I-I
-
K
+E I
~[~t[T"(,l,)
-
T"(r/,)]
T~I~(,1,)
(26)
1
l:~(r)Se(r)[hQr, r x) - h " ( r , rx)]es[l(r)rdr
Parameter ident(/ication in radia/ ffou. R. 71 Hudspeth et al. "R
- t
-
-
[drlS~'(rIS"lrwdr [~lr) T~lr)T~',lr)
r"
(T"lr)= (T~"lr)]rdr (r
+
i
,:.~',[Hh(t}-- hPlr,tl]wlr, Hdt
where JcjIp) = ~ J { p ) ~ { ")is the partial derivative of(18) with respect to the parameters :~ and fl in {231 and 2" in 125d). Note that the only head derivatives which are to be computed are those of the predicted head values hP(r, t) which are smooth. Derivatives of the measured head values (the sensitivity coel~cients) have been completely eliminated from the gradient equations. Partial derivatives with respect to the parameters to be identified are required to minimize the objective function {18). In both of the gradient methods {steepest descent and Fletcher-Reeves), a one-dimensional line search procedure was employed for each iteration to determine a local minimum of the objective function for each parameter group along their respective directional vectors. A global optimum is defined as a set of all of the parameters which minimize the objective function with respect to the minimization criteria.
4.2 Steepe,vt ek',s'cent gradient method The steepest descent method is one of the oldest and simplest gradient methods used (Bazarra and Shetty~). The optimum step size is a function of the parameter being identified. Selecting a single step size for all the parameters results in less computer time per iteration. Unfortunately, this often makes convergence to a global nfinimum dependent on the initial guess. If a different step size is used for each parameter, dependence on the initial guess is reduced. Once optimum step sizes have been estimated, new estimates for the coefficients of the parameters in (23)and the leaky parameter in 126) may be computed from
~alue of 126) for the kth iteration. When iteration k is equal to an integer multiple of the total number of parameters being identified, the steepest descent method is applied. This pacer method insures that the parameter identification algorithm does not diverge from the global minimum (Bazarra and Shetty ~1. After the gradient vectors (29) are determined, the line search technique is again applied to locate the optimum step size for each parameter group as before (Bazarra and Shett,v ~ and Gill, et al.~e). Our results confirm the conclusions given by' Neuman :~ and Carrera and Neuman 5 that the conjugate gradient method of Fletcher and Reeves improves the convergence rate of the parameter identification algorithm. Searching for a global minimum continues until the following criteria are satisfied: 1. The maximum number of iterations specified by the user is reached. 2. jk(p) <~ 10 s 3. iAJkip) J k ~(p):! ~ I0 s for 4 consecutive iterations. 5. M O D E L V E R I F I C A T I O N U S I N G S Y N T H E T I C DATA Synthetic data with known parameter distributions with and without random noise were used to test the algorithm. Synthetic data tests the ability of a parameter identification algorithm to identify a known distributed parameter system (Carrera and NeumanS). These synthetic tests may be grouped into the following categories: (1) (2) {3) (4) (5)
Identification of analytical parameter distributions; Stability of the parameter identification algorithm; Uniqueness of the parameter estimates; Effects of scaling on the identification algorithm; Comparison between the steepest descent and Fletcher-Reeves gradient methods; (6) Effects of noise on the identification algorithm.
5.1 hlent(fication of analytical parameter distributions "Measured' piezometric heads h"(r, t) for the two cases
,'= Ii! +::!)I!i where ~.~ is the optimum step size and d~.~is the directional vector for the kth iteration for parameter (-)[Bazarra and Shetty~).
4.3 Fletcher-Reeves gradient method The Fletcher-Reeves rnethod is based on conjugate gradients for a quadratic equation (Bazarra and Shetty'). The steepest descent method generates the first iterate value. For subsequent iterates, a positive multiple of the search direction G~.j from the previous k - 1 iteration is applied to the directional vector for the kth iteration according to
analyzed were computed numerically by the transformed BVP (17). Dimensions and analytical parameter distributions used for the two cases are listed in Table 1. The number of coef'ficients used to estimate the transmissivity for Case I listed in Table 1 by (23a) was N, = 2 and storativity was specified. The number of coefficients used to estimate transmissivity for Case2, Table 1 by (23a) was .,\'~ = 4 and storativity by (23b) was N, = 2. The following Root Mean Square (-RMS) errors ac~ were used to quantify' errors for hlr, t), T{r), and S(r): .
. 1 . x . 1 ~t
~ = ~
m ~] -- hi..)-
[4o
n=l
,
[(h~
L ~
(30a)
i=l
,L [ITI'-TT")-']
a~-= 7 2-
"F-
t30b)
i=I
where
, 'I k
G~.) - -
I2
, -~[ ~~l.~ - - 1 1 ,, 72
I ,L [ ( s [ - s D - ' ]
j~.~ .. o~: . . -.
I.I
(29a,b)
where I(')iJ is a Euclidean norm o f ( ' ) and J~.~ is the
(30c)
i=I
where F/° is the average value of the initial head data H°{r) in (6a); T and S are the average values of Tim and S'~", respectively: I = 25 is the total number of simulated
Adv. Water Resources, 1991. I'~d. 14, No. 5
245
P~lrameter Me, t(1ication in mdial f l o w R. T. Hudspeth et al. 30
25,
25
20
2O
]
I
2" 5
0
~0 F.--
5¸
J
" o-~
0
(a)
0
0
2 4 6 R~zdiaI Oistance
8
1o
0
(a)
20
0.5
10 Rodial
I5 2_.0 Distance
2.5
20
15 >,
~5
!
to ~n
c t3
21,22,25,24,25
IO
E c ID
I0 Theoretical
Theoretical Values
5
(b)
l
o
I
z
i
I
,~
I
i
1
6
i
8
5 ~
7
0
0.5
Vol
,
28,29,50
i
Io
Rodic]l Distance
Fig. 2 Transmissivity estimates for Case I, Table 1 with storativitv speciJied. (a) k = 0-4, (b) k = 21-25 I
measured heads; N = 21 is the total number of discrete time values. Results are illustrated in Figs 2 through 4. The solid triangles on these figures represent the location of the "measured' prior data for T(r) only for Fig. 2; and for both T(r) and S(r) for Figs 3 and 4. In these cases and all others examined, the parameters were well-estimated by the identification algorithm. In most cases the maximum number of iterations criterion # 1 . in § 4.3 was reached before either one of the other two criteria in § 4.3. These results did not depend on whether the simple basis functions (23) or more complex representations were used to estimate the parameters. 5.2 Stahilitv o/ the parameter identification algorithm Stability is achieved, in part, through the adjoint BVP (25) by eliminating the need to differentiate the heads with respect to the parameters. The line search procedure
246
Adv. Water Resources, 1991. Vol. 14, No. 5
(b)
t
I
I
1.0 1.5 2_.0 Redial Distance
2_.5
Fig. 3 Tra,smissivity estimates Jar Case 2. Table 1 • (a) k = O~l. (b) k = 21-25 and variable step size were utilized until one of the convergence criteria in § 4.3 was satisfied. Since the analytical values for the distributed parameters listed in Table 1 were given, convergence was usually obtained during the first four iterations. The program would normally terminate at this point having found the global minimum. To test for stability, however, the algorithm was forced to continue for a total of 50 iterations with the line search procedure omitted and with a step size ~k ~(.~ = 0.1 in (27). Figure 5 compares the initial iterations and final five iterations for transmissivity only for Case 1, Table 1 with storativity specified. The distributed parameter estimates T~(r) in (23a) oscillate about the true analytical value for
Pmametcr lljdllg!/i('d[[¢~ll in radia/jlo~ l R. TI Hu~gpeth et al. 20
J
40
15
>. 30 F. 2
I
>
I, 2,3,4
20
(.,q
c
eorelica[
Values
I0
al Values o.,
0
(a)
I
t
0.5
0
1.0
Radial
z
i
1.5
2.0
2_.5
Distance
50
4o
(a)
0 0
1
i
2
Theoretical Values 30 >.
,
•
t
~
i
4 6 8 Rodiol Distance
I0
£0
27, 28,29,3C
L; *6 20
2 03
15 I0 E c O
l
o (b)
0.5
l
,.o
Radial
i
t
1.5
2.o
a
e.5
~
I0
I..-
"-- 4 6 , 4 8 , 50
,Theoretical Values
Ois~'ance
F(g. 4 Storativitv estimates f o r Case 2, Tahle /." (a) k =0-4. /h) k = 2 / - 2 5
Case 1. The magnitude of these oscillations does not increase as the number of iterations increases proving that the identification algorithm is stable. The magnitude ofthese oscillations is a function of the step size ~.~ chosen in (27).
5.3 Uniquenes.v of parameter estimates Introducing prior data T"(q 0 and S"(~t) contributes to unique parameter estimates. The effects of prior data for only T"(qt) with S(r) specified were examined by the following tests listed in Table 2: Test A and B: Case 2, Table I with S(r) specified and Dirichlet boundary conditions (2 = ~ ) and prior data for only T"bh) simulated at location q~ = 0.3, 762, and 1,524 meters. Test C: Case 2, Table 1 with S(r) specified and Dirichlet boundary conditions {,:. = :c) but only one value of prior data T"Oh) was simulated at qt = 0.3 m at the well boundary.
01
0 (b)
i
i
2
L
i
4 Radial
n
~
L
6 Distance
l
8
i
Io
Fig. 5 Stability oJ transmissivitv estimatesJbr Case 1. Table 1 with storativitv specifi:ed: (a) k = 0-4, (h) k = 46-50 Tests A and B compare the sensitivity of the identification algorithm to initial guesses for parameter estimates. Test C illustrates the need to have prior data T"(th) specified away from the well boundary. Storativity S(r) was specified for all three tests listed in Table 2. The identification algorithm closely approximated the analytical parameter distribution for Case 2, Table 1 for T(r) in both Test A and B independent of the initial guesses. However, the identification algorithm failed to approximate the analytical distribution listed for Case 2,
Adv. Water Reso,trces. 1991, Vo/. 14, No. 5
247
Parameter ideptti/ication in radial flow. R. T HudYpeth et al. TABLE 1. Summary of Test Cases Analyzed CASE i
I 1
i
k
/~S
/2/°
(m)
(m)
(m)
(m~
(m~/d)
0.I52
15,24
152
270
152g
:"
/-
Ss
(:o'2
(m~Id)
~.0 4o+1oo(
2
I
0.152
1524
609
270
(m~/d)
2.0
1.0
P/R
)
2.0
2196
8.0
40+ exp{4.6ff/R)}
I+9( ://~ )
TABLE 2. Summaryof Uniqueness Tests for T(r) only with S(r) Specified (CASE 2, TABLE I) prior da~ ]ocadon:
TEST
l ]
Itc~tion k = 0
Iteration k = 30
(m)
i
i
~,
J(P)
,~,
I
.3
1
I0
I
5.14
2
762
2
0
2
-.97
A
J(p)
361 3
B
C
1524
3
0
3
.05
4
0
4
1.75
1
0.3
l
1.0
2
762
2
0
3
1524
3
0.3
I
RMS ERROR
~l
1
5.06
2
-.57
0
3
.17
4
0
4
1.58
5622
I
1.0
1
5.03
2
0
2
-1.99
3
0
3
-.27
4
0
4
-.17
5476
~k
qr
(~)
(~)
0.61
.2
l.g
.63
,0~
1.7
.40
.I3
27.9
TABLE 3. Effects of Scaling on Objective Function Terms (Case I, Table 1) TEST
(m)
(mild)
A
1.0
1.0
[
1.0
(m)
1.0
(d)
1.0
B
270.
40.
I
I0 "i
1524
2.0
]tcr'zdon k
-:(p)
J,(p)
I
I
I
I
J~)
B
A
B
A
B
A
B
o
5XlO J
3.09
5XI~
2XIO "~
5XIO 1
3
23
-0
A
1
4X10~
0.6
4XI~
6X10 "~
265
.62
.16
2
472
-0
4-45
-0
27
-0
2XIO 2
3
.27
.25
.01
1XIO 2
4
3XIO ~
2XlO"
lXl04
2XI0'
5
18XI0 "s
2XI0 ~
6XI0 7
2XI0 "~
6
16XIO s
2XIOa
6XI0 "7
--0
7
25XI0 -j
3XIO ~
6Xlff ~
2XI0 ~
6XI0 "~
8
18XIO"#
9
18XIOq
10
10XI0 "I
I"
6XlO a
2XI0 ~ LXI0 ~
Table 1 for T(r) in Test C when only one value for T"(r/t) was provided near the well b o u n d a r y (qt = 0.3 m). The parameter distribution for T(r) for Case 2, Table 1, estimated by Test C, Table 2 is illustrated in Fig. 6. N o t e that the final value of the objective function j3O(p) in Table 2 was nearly the same for all three tests; but the RMS errors for a T were not. Test C clearly demonstrates
248
:,(p)
Adv. Water Resources, 1991, Vol. 14, No. 5
-
6XI0 a
the inability of the identification algorithm to find unique parameter estimates without prior data T'(qt) being specified away from the boundaries when Dirichlet b o u n d a r y conditions (). = z ) are used. Again, storativitv S(r) was specified for all three cases listed in Table 2.
5.4 Effects of scaling on the identification algorithm
Parameter ident~ication in radial flo~ 20
io,ooo
-
R. T. Hu~£peth et aI.
-
I,O00
15 =x
g
>
E
er-Reeves
t~ c
I0
cal c o
Ioo
Z Theorelicol
Values
on
>
Io
u
',7 0
I
26,27, 2 8 , 2 9 , 3 ~ ~
0
o
(a)
I
i
0.5
I.o
Radial
0.I
i
2.0
1.5
Distance
i__
o
Io
20
i
30
L.__..
40
50
r terotions
Fig. 7 Comparison between steel)est Fletcher- Ree yes ,gradient methods
descent
and
20 ten iterations listed in Table 3. With the scaling used in Test B, Table 3, the error in head values Jh(P) was always insignificant compared to the prior data term Ja(p). For Test B. the indentification algorithm converged in just two iterations. However, the head errors Jh(P) were never greater than 10 -s because of scaling and contributed essentially nothing to the objective function (18). Identification of the transmissivity parameter for Test B was entirely the result of the prior data term Ja(P) because of scaling.
15
E I0 • Theoretical Value
I--
5--
5.5 Comparison between the steepest descent and the Fletcher- Reeves gradient methods
4../-
o-N i
0
(b)
0.5 Radial
I
1.0
1
1.5
,
2.0
Distance
Fig. 6 Uniqueness Jar Case2, Table 1 b3' Test C, Table 2. (a) k = 0--4, (h) k = 26-30
The effect of scaling on the rate of convergence may be significant. Two related scaling problems were evaluated: (1) the effect of improper scaling; and (2) the effect of scaling the basis functions. Data for Case 1, Table 1 with S(r) specified and Dirichlet boundary conditions (/. = m) were simulated for the two tests listed in Table 3. Test A has the effect of not scaling the algorithm; while Test B scales h~.,, S~, R~ and t, by their largest values and T~ by the reference value at the well. These two tests represent two extremes which could be applied. For both tests, prior data values for only T'(q~) were again specified at locations q~ = 0.3, 762, and 1,524 meters. The effects of not scaling terms in the objective function (18) are illustrated by Test A, Table 3 where the error in head values Jh(P) is two orders of magnitude larger than the prior data term Ja(P) for the first iteration k = 0. This difference in magnitude remained constant for all of the
Using Case 1, Table 1 with S(r) specified and Dirichlet boundary conditions (2 = :c) as input, the identification algorithm was run twice using the steepest descent method first and the Fletcher-Reeves gradient method second. The results of this comparison are illustrated in Fig. 7. The steepest descent method converged to a global optimum fastest.
5.6 Effects of noise on the identiJicatio,z algorithm Effects in measured heads hi.",, were simulated by introducing Gaussian white noise into the synthetic head data only. The identification algorithm searched only for the transmissivity distribution T(r) for Case 1, Table 1 with S(r) specified. Prior data values for transmissivity T,"~ were simulated without noise at locations rh = 0.3, 762, and 1524 meters. To evaluate whether the noise was truly random, the RMS error and maximum error in the 'measured' heads was evaluated using the theoretical values for Case 1, Table 1. The five tests listed in Table 4 simulate increasing amounts of noise. If the head noise is Gaussian, then the resulting RMS error a~ listed in Table 4 should be approximately 50 % of the maximum noise. For examples listed in Table 4, these conditions were satisfied. Table 4 indicates for Case 1, Table 1 with :~, = 5.0 and :~_, = 2.5 with head noise that the algorithm could identify the transmissivity as large as + 0.152 meters (Tests A-C).
Adv. Water Resources, 1991, Vol. 14, No. 5
249
P a r a m u w r ,dumiHuation in r¢Mial f l o w . R. T. HudYputh et al. TABLE 4. Effects of Noisy Head Data on T(r) Specified (Case I, Table i) M rasu rea:l Hem.ds
Itermtaon k = 0
~loise
Teat
0"n)
(%)
d (%)
t
Me, x E r r o r
i
~,.I
[
i%)
crr
o,
t
Ii
A
5.0xl0 ~
±0.003
6.0x10 a
-2.@,lxl0
6.88x10 a
.5.53x10 ~
0.016
a..99 2'.',9
3.09x.tO "~
15
B
5,Ox10 3
±0.030
5.3x10 ~
-2.6ax10 2
6.74x.I 0 ~
4.8ax.I 0"~
0,O65
4.99 2.50
0.289 t
20
C
2.Yx10 "i
+0.152
0.02~,
41.135
0.0284
4). 139
0.660
500 2 45
6.995
25
D
0.05
±0.305
0.056
0.269
0.136
1.02'7
1.058
5.12 2.50
30.a2
25
E
0.28
± 1520
0.312
1.50'0
0.,196
4.06~
2.730
aL .65
1055.8
25
,I
2.60
I
F o r Test E with head noise ___1.52 meters, the a v e r a g e p e r c e n t e r r o r in the t r a n s m i s s i v i t y e s t i m a t e s after 25 i t e r a t i o n s was only %. = 3 % .
REFERENCES I
5.7 Com'ht,vion,~
2
1. T h e p a r a m e t e r i d e n t i f i c a t i o n a l g o r i t h m successfully a p p r o x i m a t e d the t h e o r e t i c a l p a r a m e t e r d i s t r i b u t i o n s .
3
2. T h e p a r a m e t e r i d e n t i f i c a t i o n a l g o r i t h m was stable with small, b o u n d e d o s c i l l a t i o n s a b o u t the true t h e o r e t i cal solution. Stability was o b t a i n e d by the use of the a d j o i n t BVP. 3. T h e i m p o r t a n c e of p r i o r d a t a was illustrated in Fig. 6 w h e r e the correct p a r a m e t e r d i s t r i b u t i o n was not identified e v e n t h o u g h the o b j e c t i v e f u n c t i o n was m i n i m i z e d . F o r this case, prior d a t a was r e q u i r e d in o r d e r to o b t a i n a uniqt, e solution.
4 5 6 7
4. T h e effects of scaling on the o b j e c t i v e f u n c t i o n ilhtstrated that i m p r o p e r scaling results in s o m e t e r m s h a v i n g no c o n t r i b u t i o n s to the o b j e c t i v e function. P r o p e r scaling insures that all t e r m s in the o b j e c t i v e f u n c t i o n are comparable.
9
5. B o t h the steepest d e s c e n t and F l e t c h e r - R e e v e s gradient m e t h o d s c o n v e r g e d to a g l o b a l m i n i m u m .
10
6. T h e results for noisy head d a t a illustrated that the p a r a m e t e r identification a l g o r i t h m is stable and will c o n v e r g e to the t h e o r e t i c a l p a r a m e t e r d i s t r i b u t i o n e v e n with r a n d o m errors in the m e a s u r e d heads. T h e m a g n i t u d e of the o b j e c t i v e f u n c t i o n was n o t a g o o d m e a s u r e of the g o o d n e s s of fit for the p a r a m e t e r s . T h i s w o u l d suggest that several trials are r e q u i r e d for a p p l i c a t i o n s to field d a t a w h e r e the d i m e n s i o n a l i t y of the basis f u n c t i o n s and the radial scale m u s t be a d j u s t e d in o r d e r to d e t e r m i n e an o p t i m u m result.
ACKNOWLEDGEMENTS W e gratefully a c k n o w l e d g e the financial s u p p o r t p r o vided by the D e p a r t m e n t of E n e r g y ( D O E ) u n d e r the Basalt W a s t e I s o l a t i o n P r o j e c t ( B W I P ) to the N o r t h w e s t C o l l e g e and U n i v e r s i t y A s s o c i a t i o n for Science ( N O R C U S ) u n d e r c o n t r a c t D E - A M 0 6 - 7 6 - R L D - 2 2 2 J ; a n d the Office of N a v a l R e s e a r c h ( O N R ) u n d e r the U n i v e r s i t y R e s e a r c h Initiative ( U R I ) C o n t r a c t N 0 0 0 1 4 - 8 4 - C - 0 3 0 3 . W e h a v e benefitted f r o m n u m e r o u s s t i m u l a t i n g discussions with Drs. S t e v e B a k e r a n d Allen Liu.
250
M a x Err,or (m)
Adv. Water Resources. 1991, Vol. I4. No. 5
8
I1 12 13 14 15 16 17 18 19 20
Bazarra, N1. {5. and Shctt}, C. M., :\mllinuur Prowammin~l, Theory. and -II¢torirhm~. Wile}}} and Sons Inc., New York, 251 33(I. 1979 Carnahan, B., Luther, H. A. and Wilkes, J. O.. Applied :V,mwrical 3,1ethod,~, John Wile}}}' and Sons, New York, 113 115 and 429 461, 1965 Carrcra. J. and Ncuman, S. P, Estimation of aquifer paran'~etcrs under transient and steady state conditions: I. maximum likelihood method incorporating prior infc,rmation, II'RR, 22. 199 210, 1986 Carrera, J. and Ncuman, S. P. [:.stimation of aquifer parameters under transient and steady state conditions: 2. uniqueness. stability, and solution algorithms. WRR, 22, 211-227. 1986 Carrera. J.and Neuman. S. P. Estimation of aquifer parameters under transient and stead,, state conditions: 3, application of synthetic and lield data, WRR, 22, 228 242. 1986 Chavcnt, O. History matching by use of optimal theor,',, S. 0! Pet. Knelt.. 74 86. 1971 Cochran. J. C., Wiser. W. C. and Rice, B. J. ,-tdl'anced Enqineering Muthemuti{s, Brooks,Cole Publishing, Monterey. CA, 175-184, 1987 Cooley, R. l_. Incorporation of prior information on parameters into nonlinear regression groundx~ater tlov. models: 1. theor,,,. WRR, 18, 965 976, 1982 Ehlig, C. and Halepaska, J. C., A numerical study of confinedunconlined aquifers including etlkzcts of delayed yield and leakage. WRR. 12. 1175 1183, 1976 Friedman. B. Principles and Techniques Ol Applied Mathematics, Wile,,' and Sons, New York, 42 147, 1956 Gerlach, J. and Guenther. R. Remarks on part, meter identification, Convergence of the algorithm, Numerical .%lathemutics, 51, 3 9, 1987 Gill, P. E., Murray, W. and Wright. M. H. Practical Opfimizution, Academic Press. 1986 Guenther, R. B., Hudspeth, R. T.. McDougal, W., Gerlach, J. Remarks on parameter identification, I. The algorithm. Numer. 3,lath. 47. 355-361, 1985 Guenther, R. B. and Lee, J. W. Partial D!fferemial Equations of .~,lathematical Physics and Intewal Equations, Prentice-Hall, Englc;vood Cliffs. New Jersey. 16, 19 and 362, 1988 Hadamard, J., Le Prohl&ne de Cauehy et les Equations aux Dt;rivdes Partielh's Lineaires ffyperholiques, Hermann, Paris. 1932 Hildebrand. F. B. Finite-Dffl~,rence Equations and Sim,dutions. Prentice-Hall, Inc., p. 221. 1968 John, F. Partial differential equations, Ch, V, Mathematics Applied to Physics. Springer-Verlag. Ne,a York, NY, 273 277, 1970 Kravaris, C. and Seinfeld, J. Identification of parameters in distributed parameter systems by' regularization. SIAM J. Control and Optimization, 23,217 291, 1985 McWhorter, D. and Sunada. D. K. Ground~tater ttydrology and Hydraulics. Water Resources Pub.. Littleton, Co. 177-194, 1977 Neuman, S. P. Calibration ofdistributed parameter groundwater tto~v models viewed as multiple objective decision process under
Parameter identification in radial J t o w
21
22
23
24
25 26 27
uncertaint), WRR, 9, 1006-1021, 1973 Neuman, S P. A statistical approach to thein,,erse problem of aquifer h?drolog?: 3. improved solution method and added perspecti,,e, WRR. 16. 331-346, 1980 Neuman. S. P. and Yakowitz, S. A statistical approach to the in;erse problem of aquifer h.,,drology: 1. theor,,, ~IRR, 15, 845 860. 1979 Neuman. S. P., Fogg, G. E. and Jacobson. E. A. A statistical approach to the inverse problem of aquifer hydrology: 2. case stud,.. WRR, 16, 33-58, 1980 Role). K. L. An indirect parameter identification algorithm in radial coordinates for a porous medium, MS 77wsis. Dept. of Cil'il Enqmeerin~l, Oregon State University, Corvallis, O R. 1992 Ruston. K. R. and Redshaw, S. C. Seepa~le aml Grolmdwater Flm~, Wile~, Nev, York, NY. 32 35. 1979 Stakgold, l, Green's Functions and Boundary V~due Problems, John Wiley & Sons. New York, NY, 43, 1979 Sun. N. and Yeh, W. W. Identification of parameter structure m groundwater inverse problem, WRR, 21,869-883. 1985
28
29
30 31 32 33 34
R. 7-. Hudspeth et al.
Theis, C. V., The relation betv,een the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage, Tr~ms. .4mer Geophys. Union. 2, 519-524, 1935 Tikhonov, A. N. Solution of incorrectly formulated problems and the regularization method, Sot. Math. Ookl., 4, 1035-1038, 1963 Tikhonov, A. N. Regularization of incorrectly posed problems, Sot'. Math. Dokl., 4, 1624-1627, I963 Vemuri, V. and Karplus, W. J. Digital Computer Treatment of Partial D(~'erential Equations, Prentice Hall. Englewood Cliffs, NJ, 281-282, 981 Willis, R. and Yeh, W-G. Groundwater Systems Plannin9 und Munu~lement, Prentice Hall Publishing, 1987 Yakowitz, S. and Duckstein, L. Instability in aquifer identification: theory and case studies, WRR. 16, 1045-1064, 1980 Yeh, W. W-G. Review of parameter identification procedures in groundwater hydrolog,', : the inverse problem, WRR, 22, 95-108, 1986
Adv. Water Resources, 1991, Vol. 14, No. 5
251