Adaptive methods for parameter identification in ground water hydrology K . - H . H o f f m a n n , P. K n a b n e r a n d W . S e i f e r t
Lehrstuhl fiir Angewandte Mathematik I, lnstitut ]'fir Mathematik der Universitdt Augsburg, UniversiMtsstra2, D-8900 Augsburg
The paper considers the simultaneous identification of parameters and source term appearing in elliptic Partial Differential Equations from knowledge of the reference potential, which satisfies Dirichlet boundary conditions. Identifying the unknowns means complying with the equation error criterion. The adaptive method consists of transforming both the potential and the unknown parameters into time-dependent quantities, which evolve according to given laws. In the continuum case, the time decay property of a Lyapunov function is shown to hold. Next, the finite-dimensional (Galerkin) approximation of the evolution equations is considered and, given the existence of a solution, the pertaining convergence results are stated s. t., the time-dependent potential approaches the reference potential and the unknowns asymptotically comply with the time-independent equation, Then, the discrete identification problem is stated. The finite element method is used to obtain the discrete equation error criterion. A Eulerian numerical scheme for all discretized evolution equations is provided. The sequences of the corresponding solutions i. e., the potential and the unknowns, are shown to converge. The detailed proof is provided. Finally, the above method is applied to the identification of hydraulic conductivity, leakage parameter and source term of a two-layered aquifer in the Aachtal (Germany) area. The numerical results obtained from noiseless and noisy potential data are provided.
1. I N T R O D U C T I O N
In this paper we outline the adaptive methods for parameter identification in distributed systems and discuss their application to problems in groundwater hydrology. It is a general method which can be applied to the identification of an arbitrary selection of all appearing parameters in a steady or unsteady equation or a system of such equations. Here we restrict to the steady situation. Finally we will apply the method to a system, modeling an aquifer consisting of two pervious strata, but for the sake of notational simplicity we will develop the method for one equation. From the field of groundwater hydrology this covers the conservation equation of a steady inhomogeneous isotropic aquifer
-v.
in
(
h=h
220
vh) +
- w) = q
f~ C N2, on
F = 0ft,
Adv. Water Resources. 1991, Vol. 14. No. 5
(i)
where h, T, w, A, q denote the piezometric head, the transmissivity, the head of some surface water, the leakage parameter for the exchange groundwater-surface water, and the source/sink term, respectively. We can allow both for a confined and an unconfined aquifer, where T = I C ( h - b)
and K is the hydraulic conductivity and b the aquifer elevation. Ilereby T (or If in the unconfined case), A and q are functions of space and all or some of them have to be estimated from the knowledge of the head h. As our method may be classified as a method aiming at an equation error criterion, it requires full information on h. The problems connected to this fact due to the sparsity of measurements have already been discussed intensively (cf. e.g. Yeh el.al, t3) and have to be overcome by some interpolation technique like kriging. We assume this to be done and will not discuss this problem furthermore. In Section 2-4 we deal with an elliptic boundary
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
value problem similar to (1). In Section 2 we develop the adaptive method for a differential equation model, whereas in Section 3 we discuss its discretization and the interpretation of the iterative method, which is achieved. Section 4 contains a rigorous proof of the discrete convergence result and may be skipped by readers only interested in the definition of the method. In Section 5 we indicate the application of the adaptive method to a two-layer aquifer, for which field data were available.
parabolic partial differential equation
O,~(z,t)
V~ . (a(z,t)V~u(z,t))
-
+b(z, t),,(z, t) = f(=, t)
(6)
for (z, t) E a x [0, oo), u(z,t)
= g'(z)
(r)
and for (z, t) E r × [0, ~ )
~,(=, o) = uo(=) 2. A D A P T I V E M E T H O D S F O R T H E I D E N TIFICATION PROBLEM IN PARTIAL DIFFERENTIAL EQUATION MODELS In this Section we present the idea of adaptive methods for determining unknown parameters and functions in a partial differential equation when the solution is known e.g. by measurements. We restrict our investigations to the steady state case which means that as a model problem only elliptic partial differential equations are considered. The whole theory is developed in function space approach. We now formulate the basic problem. To this end, let f2 C IR'~ denote some open and bounded domain with sufficiently smooth boundary F := 0f2 (locally Lipschitz). The problem we are going to study is governed by the following linear elliptic partial differential equation of second order:
-v.
(a(z)w,'(=))
+ b(=)u'(z) =
f(z)
(2)
for z E f 2 and
~,'(z) = f(z)
(3)
for z E F. The notation " means that the corresponding quantity is known whereas the functions a(x), b(x) and f(z) are unknowns. In order to formulate our problem more precisely we use the weak formulation setting. Let Hi(f2) be the standard Sobolev space and define v
:= {~ • H~(~) I ~(z) = g'(z)
The function g'(z) Hi/-~(F) and u ° • V. Find a* • L°°(f2) stant c, b* • L~(f2) such that
for
~ • r}
(4)
is assumed to be in the space The problem then reads: with a'(z) >__c > 0 for some conwith b'(x) > 0 and f" • L~-(f2)
for x E f2 and some initial value function uo(x). The corresponding weak formulation reads:
+
fa
a(x,t)Vru(z,t)TVr(v(z)
+
fa
b(z,t)u(z,t)(v(z)-u(x,t))dz
=
fa
f(x,t)(v(z)-u(z,t))dx
We have, taking v := u* in (9):
dr(t)
=
0 fa (~u(=,t)) (u(=,t)-- u'(=)) dz
= - fa a*(z)V':u(z,t) T v , (u(z,t) - u'(x)) dx -
ffl b*(z)u(x,t)(u(z,t)- u'(x))dx
+ fa f ' ( x ) ( u ( z , t ) -
u*(z)) dz
= - fa a'(z)V=u*(x) TVr ( u ( z , t ) - u'(z)) dx -
fa b'(z)u'(z) ( u ( x , t ) - u ' ( z ) ) d z
+ fa f ' ( z ) ( u ( z , t )
- u'(~:)) dz
- fa a'(z)V,: ( u ( z , t ) - u'(z)) T
= -
(~)
holds for all v E V. To attack this problem we embed it into a time dependent problem and consider the limit case when time t tends to infinity. Instead of (2) let us study the
(9)
1
-
= [ f ' ( z ) ( v - u')(=) dz Jr~
u(z,t))dx
v(t) := ~ f . (~,(=, t) - ~ ' ( z ) ) : d~.
v ~ (~(=,t) - ~'(~))
+ jo b'(z),,'(z)(,, - ,,')(z) dz
-
for all v E V, t > 0 and u(x,0) = uo(z), uo E V. For motivation of later steps let us replace in (9) a(z,t), b(z,t), f(x,t) by a*(z), b'(z), f ' ( z ) respectively and consider the Lyapunov functional V, defined by
/ n a ' ( z ) V u ' ( x ) T V ( v - tt")(Z) dz I
(8)
f.
dz
b ' ( z ) ( u ( z , t ) - u ' ( z ) ) ( u ( z , t ) - u'(z)) dz
fa a'(z) T= (u(=,t)
- u'(z)) r
v . (~,(z,t) - u'(=)) d= - fn b°(x)(u(z,
t) - u'(z)) =d= <_50
(10)
This gives us, using a'(z) > e > 0 and Poincare's inequality:
Adv. Water Resources. 1991. Vol. 14. ,Vo. 5 221
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
d 0 < v(t) < - c ~ v ( t )
- ffl b ° ( x ) ( u ( x , t )
(11)
- u'(z))
- f~ b ' ( x ) . ' ( = ) ( . ( = , t ) for some positive constant C. From (10) and (11) follows:
u(t)--*u"
in
L2(fl)
for
+ f a / ' ( = ) (.(=,t) - .'(~)) d=
t---*z~.
=
1 f n ( u ( = , t ) _ ~'(=))"d=
+
-~ i f a (~(=,t) - a'(~)) ~ d=
+
~1 fn (b(=,t)
+
-~
d= b" (~))~"
-
(f(z,t) - f ' ( z ) ) 2 d z .
- . ' ( = ) ) d=
- f~ (b(=,t) - b'(=)) ~,(x,t) (u(~,t) - u'(=)) d~
+ ff~ (f(z,t) - if(x))(u(z,t) - u'(x)) dz
Let us now return to equation (9) and study the limit behaviour of a(x,t), b(z,t) and f(z,t) for t - c~. To this end we define a corresponding Lyapunov function V(t) by v(t):
2 dx
+ ffl
-~a(z,O
t)) (a(x,t) - a'(z)) dz dx
+A + fn
Of(x,t))(f(x,t)-f'(z))dx.
(12)
It is our goal to end up with the same estimate for ~V(t) as in (10). Therefore we define the evolutionary development of a(x,t), b(x,t) and f ( x , t ) by the following set of equations in weak formulation:
Proceeding in the same way as in (10) we obtain:
=
0
fa
-~a(x,t)v(=) dx
fa
v(x)V=u(x,t)TVx(u(x,t)--u*(z)) for all
v e L°°(a),
-
ffl b(z,t)u(z,t)(u(x,t)
~tb(x, 0 4 = ) dx
f.
v(=).(=,t)(.(=,t)-.'(=))
- u*(x)) dx
+ f~ ( O b ( z , t ) ) (b(x,t)-b'(z)) dx
=-
0
-~ f ( x , t ) v ( x )
f~
v(~)(,,(~,t)-~,'(~))d~ for all
~v(o
- u'(=)) r
V . (u(x, t) - u ' ( x ) ) dz - fn ='(x)V=u'(=)Tv': (u(x,t) -- u'(z)) dz - f~ (~(=, t) - ~" (~)) v ~ ( ~ ,
t) ~
v:,, (u(:~, t) - ~'(~)) d~
Adv. Water Resources. 1991. Vol. 14, No. 5
v ~ L-~(~).
(15)
=
f~ ,r (x)v:,, (u(x, t) - ,.,.(:~)) 'r
V=(u(x,t) - u'(x)) dz -
fab'(=)(~(~,t) - ~'(~))~a~ < O.
Unfortunately the form of V(t) does not allow to obtain an estimate as in (11). But integrating with respect to t we can still guarantee that sup t>O
222
dx
conclude from (12):
( f ( x , t ) - f*(x)) dx
f . a'(=)v~ (~(=,t)
(14)
Taking (13), (14), (15) and (5) into consideration we
-
v G L~(fi),
f~
+ fa (-~a(z,t)) (a(z,t) - a°(z)) dz
=
d=
( ~ ( ~ , t ) - ~'(=)) d~
+ f~ f ( x , t ) ( u ( x , t ) - u*(x)) dx
+ f~ ( ~ f ( z , t ) )
(13)
f~
for all = - f . ~(~, t ) v . ~ ( ~ , t ) ~ v .
dx
V(t) < oo
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
holds for all vh E Vh. Then (16)-(19) have a solution uh(z,t), ah(X, t), bh(x, t), fh(x, t) satisfying
and
f~o L v~ (~(.,t) - . ( ~ ) ) r V~ (u(x,t) - u'(z)) dzdt < c~ hold. Turning now to finite dimensional Galerkinapproximations of (9), (13), (14), (15) we will be able to prove convergence as t --* ec. To this end let Vh, Ah, Bh, Fh be finite dimensional subspaces of V, L ~ (f2), L ~ ( ~ ) , L2(fl) respectively and consider the system of ordinary differential equations: 0
/n ( - ~ U h ( Z , t ) ) (Vh(X)-- Uh(X,t)) dz
+ L ~h (~:,t)V~uh(:~, t)r V~(vh (~:) --Uh(X,t)) dx
+ L bh(~,t)u~(~,t)(,~(~) - ~(~,t)) d. = ~fh(x,t)(vh(z)-
uh(z,t))dz
(16)
for all Vh EVa, t > 0 and Uh(x,O) = Uo(Z), Uo E Vh,
~h(x,t) a~(~,t) bh(~,t)
- u'(~) - ~=(=) - - b=(=)
fh(x,t)
-- f~(x)
as t ---* oo, where a m, b~ , f ~ and u" solve
(20). This theorem says that we have to solve the ordinary differential equation system (16)-(19) for t large in order to obtain a "good" finite dimensional approximation for the unknown parameters a', b', f° in (20). The assumptions of the theorem can be replaced by some special choice of the subspaces Vh, Ah, Bh, Fh. For details see Reference 11. The differential equations (16)-(19) are nonlinear with respect to the unknown functions. This disadvantage can be avoided by following a modification of the main idea. This improvement was proposed by Baumeister and Scondo 2. In their paper they substitute the system (16)-(19) by
f. ~,~(~, t)v~(x) d.
f,~ ~uh(x, t)(vh(~) - u,,(x, t))d:~
= f~ v~(x)v~u~(:~, t ) ~ v . ( ~ (~, t) - u'(x))d~
+ f~(cv.(u~(~,
(17)
t) - u'(~))
+a(x, t)V~u'(x))r V.(vh(~:) - uh(~:, t)) dx for all Vh E Ah, t > 0 and ah(x,O) = ao(x), ao E Ah,
+ f~(~(uh (x, t) - u'(x)) + bh(~, t)u'(~))(v~ (~)
f~ ° bh (x, t)vh (x) dz =
fn V h ( X ) U h ( X , t ) ( U h ( Z , t )
--
--Uh(X, t)) dx
U*(X))dx,
= f f l h ( x , t ) ( v h ( z ) - uh(x,t))dz
(is) for all Vh E B h , t > 0 and b h ( x , O ) = bo(x), bo E B h ,
L ~fh(~, t)~(~)d= = f~ v.(~)(u'(~)
-
u~(~,t)) d~,
(19) for all vh • Fh, t >__0 and fh(z,O) = fo(x), fo E Fh. The initial functions uo(z), ao(x), bo(x), fo(x) can be arbitrarily chosen. The system (16)-(19) was considered first by Alt, IIoffmann, Sprekels 1 (see also Refs 7, 8, 11). The proof of the following theorem is contained in Ref. 11. It is an extension of the ideas in Ref. 7. Theorem 2.1 Let u ° E Vh and assume that parameter functions a" E Ah, a*(x) _> c > 0, b" E B h , b'(x) >_ c > 0, and f* E Fh exist such that
fa
a'(x)Vu*(x)Tv(
vh -- U*)(X)
fa f'(x)(vh
-
u')(x)
for all vh E Vh, t >__0 and Uh(x,O) = uo(x), Uo E Vh. C is any symmetric positive definite matrix and o" an arbitrary nonnegative real number. (See Section 3 (43) and (46) for the role of C and or).
fa O ah( x' t)vh (x) dx
= ~ v,~'(~)~V.(u~(~,t) - u'(~))~h(x)d~ (22) for all vh E Ah, t >__0 and ah(x,O) = ao(z), ao E Ah,
f~ ~ b h ( ~ , t)~(x) a~ =
ffl u ' ( z ) ( u h ( x , t) - u'(z))vh(x) dx
(23)
for all vh E B h , t >_ 0 and bh(x,O) = b0(x), bo E B h ,
dx
In ~ f h ( z , t ) v h ( x ) d z
+ fa b'(z)u'(z)(Vh - u')(x) dx =
(21)
=
fd~'(~)
-
~(~,
t))vh(~) d~,
(2o)
(24)
Adv. Water Resources. 1991. ~271. 14. No. 5
223
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
for all t'h E Fh, t > 0 and fh(x, O) = fo(x), fo E Fh. (See Section 3 (42), (44) and (46) for the role of a0, b0 and f0.) Proceeding in a similar way as in the previous theorem convergence for t ~ ¢c can be proved. Due to the linearity of the equations the proof turns out even simpler.
Theorem 2.2 Let u" E Vh and assume that parameter functions a" E Ah, b" E Bh and f " E Fh exist, satisfying (20) 1. Then (21)(24) have a solution uh(z, t), ah(z, t), bh(z, t), f ( z , t) with the property
~(~,t)
L2(f2) are given, where If" is some closed and convex subset of H*(9..). Then the problem we want to treat now reads'. Find a" E K such that, for any t, E K ' ,
fn ~ ' ( ~ ) w ' ( ~ ) r v ( u " - ,(x))d~
Our method relies strongly on the assumption that (25) has a solution. Now we consider the finite dimensional approximation of problem (25). To this end let Vh C H~(f2) and Ah C L~(f2) denote finite dimensional subspaces satisfying the compatibility relations u" EVh
--u'(~)
(25)
- f~z f ' ( x ) ( u ' ( x ) - v(~:)) dx <_ O.
and
(N'NVh)-u':;k
{0}.
(26)
Furthermore let us assume that
A(z,t)
---~foo(x)
P(ic) < I,',
as t tends to infinity and aoo, b°°, f~o and u" solve (20). For numerical purposes it is rather more convenient and it shows better results to use this adaptive method to identify the parameters a', b', f" than to apply the one we have presented first. T h a t is the reason why we restrict our further investigations in the following Sections to this technique only. The adaptive methods for identifying parameters admit numerous useful extensions. In Reference 8 this method was developed for variational inequalities in the stationary as well as in the evolutionary case. As an example serves the identification of the porosity matrix in a problem describing the saturated-unsaturated flow through a dam. The method can also be modified to become applicable for nonstationary parabolic partial differential equations. It is also possible to include restrictions on the solution and on the parameters. This might be useful if some additional apriori information on the solution of the problem is known. We shall shortly sketch a result we obtained in this situation. Let K be a set of physically admissible coefficients which is assumed to be nonempty, closed and convex, and with some c > 0,
I':C{aEL~(fi)la(x)>_c>O
a.e. onfi}.
Next we assume that measurements u" E K*, f" E
(27)
where P stands for the L2(f2)-orthogonal projection map onto Ah, and that V.Uh(Z) r . V.(Uh(Z) -- Vh(Z)) e dh
(28)
holds for u, vh E t(" f'l Vh. Our finite dimensional analogon of (25) is: Find a" E K such that, for any vh E K" N Vh,
f~ a ' ( x ) v u ' ( ~ ) r v ( ~ ' ( ~ ) - ~ ( ~ ) ) d ~ -
fnf'(x)(u'(z)
- Vh(x))dz <_ O.
(29)
For this problem the adaptive method can be successfully formulated as well.
Theorem 2.3 For any (Uo, a0) E Vh x (K M Ah) there exist unique ah E Hl(O,-~o;Ah), uh E HI(o, ec; Vh) such that ah(x,O) = ao(X), Uh(Z,0) = Uo(Z), and, for a. e. t > 0,
f~ ( ~ ( x , t )
- f.(~))
(Uh(Z, t) - vh(z)) dz + faah(z,t)V.u(z,t) T V~(uh (z, t) -- Vh (x)) dx < 0
(30)
for all vh E K" M Vh,
ffl (~tah( x, t)) (an(x, t) - dh(z)) dx
- fn(ah(~, t) - d h ( = ) ) V . u h (x, t) T 1 W i t h o u t this a s s u m p t i o n we can still prove t h a t the solution of (21)-(24) becomes s t a t i o n a r y as t tends to infinity. These s t a t i o n a r y limits a ~°, b°°, foo a n d u °° t u r n out to be an approximate solution of (20) in some least-squares-sense. These limits are identical with the limits we get for the fully discrete scheme (41) in Section 3. Look at Section 3 (46) for a characterization o f a ° ° , b °o , foo a n d u ~ .
224
Adv. Water Reso,o'ces. 199/. Vo/. 14..Vo. 5
V,:(uh(x, t) - u'(Je))dx <_ 0 for all dh E [( I'1 Ah. Furthermore:
~h(x,O ah(~,t)
--.~'(~) --,aoo(x)
(31)
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
as t
--
c~
and a~°(z) solves problem (25).
Since this result is not applied in this paper we omit the proof. The details together with further results will appear elsewhere. In the next Section we turn to the discrete problem which is the basis for computational applications.
problem. The dimension t: should be chosen sufficiently large to properly approximate the boundary value problem. The discretization of the unknown functions is accomplished by looking for them in finitedimensional spaces, i.e. one chooses e.g. linear independent functions
¢° 3. A D A P T I V E M E T H O D S F O R T H E CRETE IDENTIFICATION PROBLEM
to span the space for f,
(aVu', 27¢,) + (bu', 5,) = (f, ¢,) i = 1,...,k.
(32)
Ilereby k < t: and ¢ k + 1 , . . . , ¢ ~ correspond to nodes at the boundary 092, where for the boundary value problem the values of u are fixed (the Dirichlet condition). We use the abbreviation
(1, g):= (V f, V9) :=- fa(vfTVg)(x)dx.
to span the space for b,
DIS-
In Section 2 we formulated the equation error criterion (5) for the differential equation model and discussed the adaptive methods to achieve a solution. To arrive at a numerical algorithm both the problem and the adaptive method have to be discretized. In turn we have to discretize both the functions to be identified, e.g. the function a, and the known observations u'. For the second aspect, it is natural to use the finite element method to approximate the underlying boundary value problem (2). Let ¢1, . . - , ¢~: be a basis of the finite-element-space for u. In our numerical examples in Section 5 we will use linear elements for a given triangulation of the domain fi, which is assumed to be polygonal. Then (5) is substituted by
for
to span the space for a,
(35) if the three of them are to be identified. Thus we arrive at the discrete equation error criterion Kaa + [Qb + Ix'If = 0, (36) where the matrices I{~, IQ, K! are defined by
(I(a)ij : = ( ¢ ] V u ' , V ¢ i ) , j= 1,...,g, (Iib)ij : = ( ¢ ~ u * , ¢ i ) , j = 1 , . . . , m , := -(0:,¢i), j = 1,...,p, and
i= 1,...,k. (37)
Analogously to u*, a, b and f have to be understood as the vectors of components with respect to (35). If not all components of a, b, f are to identified, then (36) has to be interpreted in the following way: The unknown components of a, b, f are collected in a vector v, the known components are collected in a vector w*. This allows for the fact that e.g. some components of a are unknown and some are known. The columns of Ka, Kb, and K! are assembled according to the ordering of v and w* into matrices K and L. E.g., if only a is to be identified, then
(33)
u* is now assumed to belong to the finite-element-space for u, i.e. In general (36) is equivalent to
u'(x) = E u~¢i(z).
(34)
i=l
Here and in the following we do not distinguish in notation between a function in a finite-element-space and the vector of its components, i.e. for Lagrangian elements the vector of its nodal values. This means that (for Lagrangian elements) the knowledge of u" is only required at the nodes, but these values are used not only to approximate u* via (34), but also to approximate Vu" via
=
i=1
For perturbed data this is a source of error amplification reflecting the ill-posedness of the identification
Kv
= y
(38)
with y = - L w ' , i.e. y is known. The discrete equation error criterion requires to find an (approximate) solution of (38). In early work on the equation error criterion the same discretization for u and a has been used, i.e. e = t:, Ca = ¢i (of. Ref. 4), but now it is fairly well understood that a low dimensional discretization of a, b, f has a regularizing effect, i.e. enhances the stability properties of (38) (cf. Ref. 12). The two popular approaches to choose this parametrization of a, b, f are the interpolation and the zonation method. In the first case a, b, f are chosen from some Lagrangian finite-element-space and then the parameters are the nodal values. Usually the corresponding grid is a coars-
Adv. Water Resources. 1991. Vol. 14. No. 5 225
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
ening of the grid underlying the discretization of u" (cf. e.g. Refs 3 and 9). In the zonation method ~ is subdivided into finitely many subdomains, the zones, on each of which the functions to be identified are assumed to be constant. Here the parameters are the values corresponding to the zones and the basis functions are the characteristic functions of the zones. In Section 5 we will use this approach. Also if the number of parameters is small, in general we can expect neither solvability of (38) nor, in case of existence, uniqueness of a solution. In an appropriate choice of dimensions the number of parameters will be less than the number of equations. Then the solvability of (38) cannot be guaranteed, not even by the solvability of the continuous identification problem (5) due to the appearance of the discretization error by the substitution of (5) by (38). In addition we can only expect noisy data u ' . Thus we have to substitute (38) by (Iqv - y)m M ( K v
- y) --+ min
continuous equation error criterion (5)
algorithm to solve (5)
adaptive .~ method
I discretization
I discretization discrete equation error criterion (38)
algorithm to solve (40)
adaptive --~ method
Fig. 1. Development of the algorithm
((~-'~.(u7+1 - u}')¢s)/At , ¢,) j=l +
-
(39)
j=l l
where the symmetric, positive definite matrix M in (39) may be chosen according to different criteria (see below). In Ref. 13 this approach has been applied assuming the unique solvability of (39) (for identification in an unsteady situation). But the solution of (39) need not be unique due to the possible linear dependence of the columns of K. Consider e.g. the following extreme case for the identification of a: If u" is constant on some subdomain (~ of f2 and the support of ¢ ] for some j, i.e. a zone in case of the zonation method, belongs to ~, then the j t h column of K~ is identically 0 according to (37). Therefore we have to extend (39) to:
+
(E
ajn+l.ar--, ~ivu . ,V¢i)
j=l
j=l m
n+l
+
=
i = 1, j='
for
u']+'=g;
j=k+l,...,L-,
((~--~(aT+' - aT )¢~])lAt, ¢~) j----1
( K v - y)T M ( K v - y) - - rain,
=
select the one for which
(vu'Tv(
% j =l
(v - vo) T N (v - vo) ~ min,
i= 1,...,g,
(40) where the symmetric positive definite matrices M, N and the parameter vector v0 are to be chosen. We will develop an iterative method for (40) by properly discretizing the adaptive method (21)-(24). Thus we will describe the adaptive method in the lower branch of the diagram in Fig.1. by discretizing the adaptive method in the upper branch. We discretize (21)-(24) with respect to the space variable by the finite element method as defined in (32) for u and analogously for a, b, and f in the space defined in (35). With respect to the time variable we use the fully implicit Euler scheme. To simplify the notation, let us take the time step size At to be constant. If all parameters have to be identified, this leads to
Adv. Water Resources. 1991. Vol. 14. No. 5
,¢,1
j:l P
From those parameter vectors v, for which
226
b "
( ( ~ ( b ~ ' + ' - b'])¢~)/At, ¢~) j='
i =
(u" ' ~
(. n-}-I -
•
b
j=l
i= 1,...,m, P ((Z(s;+, j=l
_
=.,,SV'(".+' u;)¢j,¢[), j=l
i= 1,...,p.
u °=uo,i,
i=l,...,k,
a 7 : no,i,
i = 1 , . . . , g,
k
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
bi° = b o , i , f O = fo,i,
i= 1,...,m, i= 1,...,p.
Hereby the superscript n indicates t '~ = n A t . The initial conditions in the (finite-element-)spaces for u, tively. In matrix notation (41) can
(41)
the nth time level are assumed to be a, b, and f respecbe written as
N~ ~"+'-~" + D(u '~+x - u ' ) At
The matrices N~, D and N~ are symmetric and positive definite and thus the same holds true for the symmetric part of the matrix on the left hand side of (45). Therefore this m a t r i x is nonsingular, i.e. for given u '~, v'~ there are unique u '~+1, v '*+1 solving (45). This solution can be computed by a variant of the conjugate-gradient method (cf. e.g. Ref. 6). An iterative method to compute u '~+1, v '~+1 is of advantage, as with u", v n often a good initial guess is available. For the sequence u '~, vn we shall prove the following result:
+K,,a '~+1 + Kbb '~+1 + K l f n+l = O, (i)
Xo
_
-
'-b" Nb b " +X-7
_ Kbr(u +
_ u'),
u'),
(ii)
= K]'(u"+:-u'),
u° = u o ,
a°=ao,
b° = b 0 ,
fO=fo.
(42)
The matrices are already defined in (37) and by
(N~)ij : = ( ¢ / , ¢ , ' )
i,j = l,...,g,
(Nb)ij := (¢~, ¢~),
i,j = l,...,m,
:= ( ¢ : , ¢ : ) ,
(iv)
D(u"
i,j = 1,...,p.
(43)
N~ ~"+'-~" ~, + D(u n+l - u ' ) + K v "+1 = y, Nv v"+i-v" K T ( u n+l - u ' ) , At yO ~
:= v
Kr(
If not all parameters have to be identified, we proceed as above and (42) transforms into
UO ~
v °° is the orthogonal projection of v ° onto the set of solutions of (ii), i.e. onto c := {,, I K v = 9} with respect to the scalar product defined by
i,j = 1,...,k,
(Na)ij :-- ( ¢ ; , ¢~),
U0 ~
Z = (Zl,..-,zk)T
(iii)
i,j= l,...,k
(D)ij := ( C V C / , V ¢ I ) + (o'¢j,¢i),
u nn~__~uoo for some u ~°, vn n - ~ v O o for some v°°, i.e. the adaptive system becomes stationary, K v °° = ~ and ~) is the orthogonal projection of y onto r a n g e ( K ) with respect to the scalar product defined by (y, Z)D-, :-- y T D - I z for Y = ( Y l , . . - , Yk) T,
VO"
(44) N~ is a block diagonal m a t r i x composed from the columns and rows of Na, Nb, and NI, corresponding to components to be identified. E.g., if all components of a and b have to be identified, then 0
-
u ~)
N
w,
= Kv ~
-
V,
" - u °°) = 0
Before coming to a proof, we shall discuss various aspects of (46). The result is more general than indicated by its derivation. The m a t r i x K and the right hand side y in (38) need not be given by (37), but m a y be completely general. In this way (45) defines an iterative method to select some least squares solution of K v = y in the sense of (40). The notion of solution is defined by the matrices D, N, and the vector v0. The matrices need not be defined by (43), but may be general symmetric positive definite matrices. The same holds true for the matrix Nu. In addition N~ and u0 do not affect the solution v °° (but the "artificial state" u°°). With respect to the choice of the scalar products, i.e. of D and N~ we have the following elementary observations: In an indirect method, using the output error criterion, usually one ends up with a pair u, v fulfilling exactly the discrete state equations, i.e. K(~)V = y,
In system form (44) reads
--AtK T o
Nv ,]
(..+,) v n+l
(47)
where K = K ( u ) indicates the affine-linear dependence of I£ on the state according to (37). In an equation error method this is not the case. If we set ~ := u*, := v ~ we fulfill (47) only up to the equation error Y-Y
:
K ( u ' ) v °~ = U + ( / - y,
,,"
u°=u0,
(46)
v° = v 0
(45)
where ~ - y is connected to u" - u ~ by (46)(iv). If we set g := u °°, g := v °~ we conclude from
Adv. Water Resources, 1991. Vol. 14. No. 5
227
AdaptWe methods for parameter identification in ground water hydrology: K.-[[. Hoffmann et al.
we get:
(46)(iv): ~ ; ( u ~ ) ~ ~ = y + d,
AtD(-ff - u ' ) = A t ( g - fl) = A t ( g - I('~) AtliT(-ff - u ' ) = A t K T D - I ( y - !t) = 0
where
(53)
d~ = ( ( c - ~ = ) v ( ~ " - u=), v¢~) + ( ( ¢ - w ) ( u " - u ~ , ¢~),
i = 1 , . . . , k,
(48)
if all components have to be identified, tIereby l
t+m
i=1
i=l+l
With the help of (53), equations (52) reduce to N~6u '~+l + A t D S u '~+1 + A t t ( 6 v ~+1 : N~6u n - A t K T 6u '~+1 + N~6v '~+t : Nv6v ~
(54) Therefore a possible strategy to incorporate prior information, e.g. a good guess of v ~ , is to choose C to be this guess. Another way to use prior information, also of statistical type, which is reminiscent of a corresponding output criterion approach (cf. e.g. Refs 10 and 5) is to choose e.g. D as the covariance matrix of the residuals (compare Ref. 13). In the same way prior information on the parameter can be incorporated by v0 and the choice of N,.
4. P R O O F OF T H E GENCE RESULT
DISCRETE
We multiply the equations (54) by 6u ~+l and 6v ~+I, respectively, and add them. This leads to: 5u ~+lT N~ 5u ~+1 + 5v ~+ffr N~ 5v~+1 = - A t 6 u n+lT D 6 u n+l A t 6 u n + l r I ( S v n + t + 6[tn+lT.:\rusun
_
+ A t 6 v n + I T K T 6 u ~ + I + 5V~+IT,'v~sv" = - - A t 6 u ~+IT D 6u ~+l + 5U"+ITN~Su ~
Using Young's inequality, i.e. uT N v < l u T N u q_ 1 T -2 -~v N v
CONVERwe get from
In this Section the convergence result (46) will be proved under the assumptions: Let K be an arbitrary matrix, say of p rows and u columns, i.e. K E IRu'~, g, u0, u" E lR u, v0 E IRV, D, N~ E IRu'u symmetric and positive definite,N~ E IRY'~ symmetric and positive definite. (49)
(55)
+ 5V~+'TN~6v~
(55):
½(6Un + l r Nu 6u '~+l + 6v n+lT N , 6v"+1) <
_At6u,~+a
+½(~
r
D 5u n+l
Nu ~u~ + 5v~ N~ ~ )
We define the discrete Lyapunov functional Vn: V,~ := ~(Su "T Nu 5u" + 5v " T N , 5v")
Then the assertions of (46) for the iteration scheme (45) hold true.
V n+l '< --At6u n+lT D 6u '~+l + V '~
: : u" + D - ~ ( y
-
~)
n+l
(50)
We take any Z with K [ = 9 and define
6ukT D S u k + V °
k=l
for all 6v '~ : = v ' ~ - ~
(58)
From (58) we conclude: V ~+1 < - A t E
6u '~ : = u " - - i f ,
(57)
Then (56) reads:
Proof
Let ~ be the orthogonal projection of y onto range(K) with respect to (., .)D-~ and define ~ by
(56)
nEPq
(59)
(51) This estimate implies:
We rewrite (44) by means of 6u '~ and 6v n and obtain: s u p V '~ < V ° < c ~ N~Su "+1
+
A t D 6 u '~+1 + AtD-5
+ A t K 6 v '~+1 + A t K ~ = N~Su '~ + A t D u " + A t y - A t K T 6un+ 1 _ A t l . ( T ~ + N~6v n+l = Nv6v n
-
and
6ukT D Su k < V ° < e~
(61)
k=l
AtI(Tu *
Using (50) and (Y - Y, k)~-~ : 0 for all k E range(K),
Adv. Water Resources, 1991, Vol. 14. No. 5
¢,o
AtE (52)
228
(60)
nElN
Therefore: 6u n " - •
O,
i.e.u ~ " - ~
(62)
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
This implies (46)(i) 1st part and (iv) by" defining u ~ := (note (53)).
(., . ) N . -
We have to prove v ~: = ~.
T h e bound (60) implies:
Suppose v ~ ¢ ~. T h e n define:
Un k k~C~ UOo
V~ :=5+A(9-v for some v °° and for some subsequence (n~) C IN. (63)
L~--U
+ D(-6- u*) + Kv ~¢ = y
(64)
By means of (50) this is equivalent to : c ~ ~ : ~,
(65)
i.e., (46)(ii) holds. We still have to prove v '~ " - ~,- v ~ , i.e., the whole sequence (v n) converges to v ~ :
for
AEIR.
~xEG
for all
AEIR
(69)
Therefore we can use v-x instead of ~ in (51) and get some V~ with the definition (57). T h e results for V n are also true for Vx~. T h e e s t i m a t e (58) implies t h a t V~ is m o n o t o n e decreasing: 1 n 1 n -i(u --U,u~--U)N.+-i(v --v~, v n --V~)N. <_ V2 (70) Passing to the limit in this inequality yields: ( v ~ - v~, ~oo _ v~)N. _< (u0 _ ~, ~0 _ ~ ) u .
Let v ~ and v~° be two limit points. Then:
"~-( U0 -- U'A, v0 -- Y'A)N,,
v nk k--T---~v ~
and
v"'
v~
for some subsequences (nk) C IN, (nt) C IN.
(71)
As ~5is the orthogonal projection v ° onto G with respect to (., .)N,, we have (66)
W i t h the reasoning above, v ~ , v~° fulfill (65), i.e. they are suitable choices for g. From the monotonicity of V n due to (58) we conclude:
(y0 -- Y'A, U0 -- Y'A)N,
= (y0 -- ~, V0 -- Y)N~
+(~-v~,~-v~)N.
v~, v ~ - V~)N.
< C + (~, - V~, ;, -
for some V ~ and for the whole sequence {V'~},eN
for all
This implies
(72)
C := (u ° - g , u ° - - g ) N , + ( v° -- ~, v° -- ~)N, is a constant independent of A. W i t h (72) we can rewrite (71) to:
(,,~ -
V ~ ~.-7~ V o°
1
(68)
Because of ~ E G and v ~ E G, we have
We a p p l y the first equation of (44) to this subsequence and pass to the limit. This leads to:
N,~ ~
~)
V~)~v.
AEIR
(73)
Using the definition of ga (68), we get from this inequality:
oo
~(v t - V ) T N v ( v F - v ) (A + 1)2(v ~ - 5, v ~ - V)N. <
1 oo -~)rN~(u~--a) + ~(u =
lim V n" = V ~ = lira V TM
-
~(v 2 - V ) T N ~ ( v ~ - v )
+
~1 ( u oo - - g ) T N . ( u
C + A2(v ~ - ~, v ~ - ~)N.
1
~-~)
for all (67)
A E IR,
(74)
that is (2,~ + 1 ) ( v ~ - ~, v ~ - ~ ) N . < c
Choosing := v~'
we get from (67)
i.e. all limit points of ( v , ) are identical. This implies v ~ n - ~ , v ~ . Therefore (46)(i) holds true. To prove (46)(iii), let ,5 be the orthogonal projection o f v ° onto G = {z I K x = Kv ~ } with respect to
for all
A E IR.
(75)
For v c° ~ 6 (75) cannot be true for all A E IR. Therefore v °° = ~ and thus also (46)(iii) holds true. In (41)-(46) we assumed t h a t the time step size At is constant. An inspection of the convergence proof shows that we can allow variable time step size Atn, if it is b o u n d e d from below, i.e. At,~ _> Attain for all n and some Attain > 0 ( c o m p a r e (61) and (64)). Because
Adv. Water Resources, 1991. Vol. I4. No. 5 229
Adaptive methods f o r p a r a m e t e r identification in ground w a t e r hydrology: K.-H. H o f f m a n n et al.
the solution of (21)-(24) converges only for t -- zc to its limit point and (41) is an approximation of (21)(24), we expect convergence of un, v,~ only for large rl t,~ := }-'~=x Ate. This convergence can be accelerated by a step size control of Ate, for which there is (theoretically) no upper bound. In Section 5 we apply the following heuristic strategy: If the relative change of the parameters is too small in a time step, the next time step is increased (multiplied by 10). Such a stepsize control turns out to be decisive for fast convergence.
and h, (kf2t~Vh2).u
on on
r~, r~.
(79)
Ilere to is the thickness of the (confined) aquifer. The sink and source term q2 includes (distributed) inflow of water through Karstic formations and the withdrawal of water by pumpage. Analogous to (21)-(24) the adaptive system for (76)-(79) reads:
O -hV . ( Cl1 V ( h l ( (t ) - h'1)t + k f)l ( t ) t ' l V h ' l )
5. NUMERICAL RESULTS FOR THE GROUNDWATER FLOW MODEL "AACHTAL"
In this Section we apply the adaptive method to a groundwater model, for which field data were available. It consists of two aquifers, an upper, phreatic aquifer, and a lower, confined aquifer, separated by a semipervious layer. The aquifers are modeled by plane domains f~l and f22 with conservation equations obtained by vertical averaging and the Dupuit approximation. The exchange between the aquifers is modeled by a distributed leakage term. Thus the equation for the steady-state groundwater flow in the upper aquifer f/1 is: -V
= ~ =q2
= ql(t)
hi(t) =
in
f~l x [0,c~)
on
rl x
( C i V ( h l ( t ) - h'l) + k f I ( t ) t ' l Vh*1) . f,
on
= ql(t)
× [0,~)
hi(0) = h °
in
O h2(t) - V . ( C ~ V ( h ~ ( t ) - h~) + k f 2 ( t ) t ~ V h ~ )
. (k f lQ V h l ) + -7-(hl - h~) n~nn~ ~
= q~
in
f~l C IR~
(76)
hl (k f t t a V h ~ ) . u
= hl = ~
on on
q2(t)
in
f22 x [0,c~)
h~(t) = h~
on
F~ x [0,c~)
=
with the boundary conditions
F[,
( C ~ V ( h ~ ( t ) - h'2) + k f 2 ( t ) t ~ V h ~ ) . u
F12.
Here k f t is the hydraulic conductivity, tl is the thickness of the saturated zone of the (phreatic) aquifer, hi is the piezometrie head of this aquifer. h2 is the piezometric head of the lower aquifer, kv is a leakage parameter describing the permeability of the semipervious layer between the two aquifers and d is the thickness of this semipervious layer. qt is a sink and source term. This term includes natural replenishment from precipitation and the withdrawal of water by pumpage. (76) is valid in a domain fll with boundary 0ill = 1" F1UFIhi is a prescribed head on F~, ql is a prescribed flow over I'12, where u is the outward unit normal vector. The flow equation in the lower aquifer f12 is similar:
on
=q2(t)
(77)
F~ x [ 0 , ~ )
h2(0) = h ° O
1
in
.
in
(f~l N122) x [0,e~) in
kv(O) = kvo 0 O-~tkfl(t) = t ' x V ( h l ( t ) - h ' ~ ) T v h [
in
f~l x [0,co) in
kA(0) = kf °
0
O--~tkf2(t) = t ; V ( h 2 ( t ) - h*2)Tvh~
in
f21 13 f2~
f21
Q2 x [0,oo) in
f~
. ( k f 2 t 2 V h ~ ) + k~-(h~ - h i ) f h n a ,
= q.~
in
f~2 C IR~
(78)
°ql(t)
= h "1 - h i ( t )
in
f21 x [0,oo)
ql(0) = qO
230
f22
k v ( t ) = -~(h 1 - h ; ) . ( ( h i ( t ) - h i ) - (h2(t) - h ; ) )
~:A(0) = ~:f °
-V
f~l
Adv. Water Resources. 1991. Vol. 14. No. 5
in
al
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
=
h;
-
in q2(O) = q~
×
in
f~
(so) Here h~, h~ are the measured pressure heads, t~, t~ are the (given) thickness of the saturated zone in the upper respective lower aquifer, d is the thickness of the semipervious layer between them. Note that there is no distinction between the confined and the unconfined part of the aquifer. Due to the equation error approach the substitution of the (measured) heads h* leads to a given thickness t~ and t~. An output error approach has to deal with nonlinear equations in case of an unconfined aquifer. The discretization discussed in Section 3 leads to the discrete adaptive method analogous to (44), which is applied here. We use linear finite elements for the heads (cf.Fig.2) and the zonation method for the pa-
rameters to be identified (cf. Figs 3,4,5). In addition, based on the same finite-element-discretization, we developed an algorithm to approximate the solution of the direct, i.e. non-linear, problems (76)-(79). In the following we use estimates, gained by personal experience, of the geological and hydrogeological data of the Aachtal, a region in southern Germany north-west of the Bodensee. This data together with the finite-element-discretization of the region were provided by Dr.-In& B. Herrling (Karlsruhe). Our intention was to test the procedure in a realistic, but controlled situation. Therefore we avoided the problem of lacking measured heads and used the parameter estimates provided to create "measured" heads by solving the boundary value problem (76)-(79) with the finite element method. We present results for the identification of k fl, kv and q2. As initial values for the heads we use the computed "measured" heads, for the parameters we took a value of the correct order of
J
Fig. 2. Finite-element-discretization
Fig. 4. Zonation of kv
Fig. 3. Zonation of k fl
Fig. 5. Zonation of q2
Adv. Water Resources. 1991, ~'ol. 14. No. 5
231
Adaptive methods for parameter identification in ground water hydrology: [(.-[[. Hoffmann eL a l
magnitude. C and (r were choosen to be the constant 1. As shown in the tables we were able to reconstruct the "true" parameters from the exact knowledge of h~ and h;. (maximum of perturbation = 0.0 [cm]). In this case the discrete identification problem has a solution (since we solved (76)-(79) with the "true" parameters to get h[ and h;, these "true" parameters are a solution of the identification problem), therefore output error and equation error tend to zero (compare (46)(iv)). This solution can be expected to be unique, since the space for the parameter is lowdimensional compared with the space for h. If the perturbations are small enough (maximum of perturbation = 0.1 [em]), the estimated (computed) parameter values are still close to the "true" parameter values. Since we cannot expect the existence of a solution of the identification problem in the case of perturbed data (we only have a solution of the corresponding least-squares-problem), output- and equation error do not tend to zero (compare (46)(iv) again). The perturbations of the data becoming bigger (1.0 resp. 10.0 [cm]), also the deviation of the estimated parameters from the true parameter increases. The parameter values corresponding to small zones differ more from the "true" values than the parameter values corresponding to large zones. This is the effect of regularisation by zonation. If the perturbations are too big, we get some negative conductivities resp. permeabilities, which have no physical meaning. In our computations we have not applied any positivity constraints to the evolution law.
CONCLUSIONS We have presented an equation error criterion method for arbitrary identification problems in steady state situations based on the idea of adaptive identification. Within the limits of a method aiming at an equation error criterion our approach seems to fulfill the requirements for an identification procedure in an optimal way: It allows for any parameterization method to stabilize the problem. This leads to a discrete equation error criterion, i.e. a set of linear equations, for which the diseretization of the adaptive methods provides an iterative procedure to select a solution in the sense of a generalized pseudo inverse. Here neither the existence of a solution nor its uniqueness is required. In case of nonexistence a least squares solution is selected. The choice of the scalar product gives the possibility to incorporate prior information on the perturbation of the observation, say the heads. If the least squares solution is not unique, in the selection criterion there is again space to incorporate an initial guess of the parameters and statistical information. Ilere one has to keep in mind that the appearance of nonexistence and nonuniqueness is enhanced by numerical computations with a finite arithmetic. One could
232
Adv. ~I~lter Resources. 199I. Vol. 14. No. 5
suggest that there are other methods to compute this generalized pseudo inverse. This is true in principle, but preliminary comparison with the truncated singular value decomposition (SVD) reveals the following: Our method seems to realize an additional regularization in the sense that its results were comparable with the SVD results for the optimal (in general unknown) truncation level. This needs further investigation. ACKNOWLEDGEMENT We gratefully acknowledge Dr.-Ing. B. IIerrling, Institute for Hydromechanics, University of Karlsruhe, for providing all data of the Aachtal area and for the permission to reproduce Figs 2-5.
REFERENCES
1. Alt, H.W., ttoffmann, K.-H., Sprekels, J., A numerical procedure to solve certain identification problems, Intern. Ser. Numer. Math., 68, 11-43, 1984 2. Baumeister, J., Scondo, W. Adaptive methods for parameter identification, Methoden und Verfahren der Mathematischen Physik, Band 34, 87-116, Frankfurt, Brosowski, B., Martensen, E., (eds.), 1987 3. DiStefano, N., Rath, A, An identification approach to subsurface hydrological systems, Water Resour. Res., 11(6), 1005-1012, 1975 4. Frind, E.O., Pinder, G.F, Galerkin solution of the inverse problem for aquifer transmissivity, Water Resour. Res., 9(5), 1397-1410, 1973 5. Gavalas, G.R., Shah, P.C., Seinfeld, J.H., Reservoir history matching by Bayesian Estimation, Soc. Pet. Eng. J.. 16(6),337-350, 1976 6. tIageman, L.A., Young, D.M., Applied Iterative Methods, Academic Press, New York, 1981 7. Iloffmann, K.-H., Sprekels, J., On the identification of coefficients o[ elliptic problems by asymptotic regularization, Numer. Funct. Anal. and Optimiz. 7, 157-177, 1984-85 8. IIoffmann, K.-tt., Sprekels, J., On the identification of parameters in general variational inequalities by asymptotic regularization, SIAM J. Math. Anal. 17, 1198-1217, 1986 9. Kunisch, K., White, L.W., Identifiability under approximation for an elliptic boundary value problem, SL4M J Control ?J Optimization, 25, 279-297 10. Neumann, S.P., Yakowitz, S., A statistical approach to the inverse problem of aquifer hydrology. 1 Theory, Water Resour. Res., 15(4), 845-860, 1979 11. Seifert, W., Identifizierun 9 yon Parametern eines Grundwasserleiters mit adaptiven numerischen Vet-
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann et ah
fahren, Diplomarbeit an der Universit~t Augsburg, 1989 12. Y'eh, W. W.-G., Yoon, Y.S., Aquifer parameter identification with optimal dimension on parametrization, Water Resour. Res., 17(3), 664-672, 1981 13. Yeh, W. W.-G, Yoon, Y.S., Lee, K.S., Aquifer parameter identification with Kriging and optimal parametrization, Water Resour. Res., 19(1), 225-233, 1983
number of iterations of the cg-algorithm to solve the linear system per time step. If this number is equal to zero, the adaptive system is stationary.
true parameter values: the parameter values used to compute h~ and h~
estimated parameter values: EXPLANATION
OF
THE
TABLES:
the parameter values estimated by the identification procedure
maximum of perturbation = k [cm]: output error: h~ and h i were perturbed by random numbers, uniformly distributed over l - k , +k]
difference between u" and u ' : Hu'* - U'HL2(n )
number o/iterations:
equation error:
number of time steps
defect in the equation (bnu • ' q~:) _ ( f . , ~3)]2)1/2
(38):
(~3[(a'~Vu',V~bj) +
number o/eg-iterations:
TABLE I Identification of "KFI" ( maximum of perturbation = number ofl number ofl timeiteraI CG-itera-I steptions I tions I size 0
0
1.00E÷03
I I I
"true" parametervalues
1. 1.80000E-03
2. 3. 4. 5. 6. 7. 8. 9. 10. 11.
3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
I I
[
0.0
[cm] ):
estimated parametervalues 1.00000E-03 1.00000E-03 1.00000E-03
I outputI error
[
] equation[ error [
0.00000E+00
6.33198E-02
1.O0000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.O0000E-03 1.00000E-03 1.00000E-03
5
19
1.00E+07 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
1.04328E-03 6.38770E-04 9.29897E-04 1.00395E-03 4.08240E-04 9.98870E-04 1.10307E-03 6.37250E-04 5.67915E-04 1.01456E-03 9.22708E-04
3.46662E+01
3.27984E-02
13
40
1.00E+10 I. 2. 3. 4. 5.
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05
1.79902E-03 3.59372E-04 7.09730E-04 2.71124E-03 6.98924E-05
2.05240E-01
2.36536E-04
Adv. Water Resources, 1991, Vol. 14, No. 5
233
Adaptzve methods for parameter identification in ground water hgdrology: K.-H. Hoffmann et al.
6 7 8 9 I0 11
I 2 7 3 1 I
10000E-03 50000E-03 20000E-05 60000E-04 80000E-03 80000E-03
1.09970E-03 2.49994E-03 7.21234E-05 3.59933E-04 1.74652E-03 1.79994E-03
18
2
I.OOE+I5
i. 2. 3. 4. 5. 6. 7. 8. 9. 10. Ii.
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
8.14081E-I0
8.01109E-13
19
0
1.00E+16
1. 2. 3. 4. 5, 6. 7, 8, 9, 10, 11,
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
8.14081E-10
8.01109E-13
TABLE 2 Identification number iterations
234
of "KFI"
(maximum of p e r t u r b a t i o n "true" parametervalues
] [ ]
=
0.1
[cm])
] outputI error ]
ofl number of[ time[ CG-itera-I stepI tions [ size
I I I
0
0
1.00E+03
1. 2 3 4 5 6 7 8 9 i0 ii
1.80000E-03 3.60000E-04 7 10000E-04 2 90000E-03 7 10000E-05 1 10000E-03 2 50000E-03 7 20000E-05 3 60000E-04 1 80000E-03 I 80000E-03
1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03 1.00000E-03
O,O0000E+O0
6,35011E-02
18
0
1.00E+15
1. 2. 3. 4. 5. 6. 7. 8. 9. I0, 11.
1 3 7 2 7 1 2 7 3 1 1
1.82166E-03 3.61354E-04 7.15133E-04 2.94227E-03 7.08022E-05 1.10933E-03 2.54685E-03 8.19355E-05 3.61599E-04
1.85093E+00
1.59893E-03
Adv. Water Resources, 1991, Vol. 14, No. 5
80000E-03 60000E-04 10000E-04 90000E-03 10000E-05 10000E-03 50000E-03 20000E-05 60000E-04 80000E-03 80000E-03
estimated parametervalues
:
1.81871E-03 1.76561E-03
I equationI error
I
Adaptive methods for parameter identification in ground water hydrology: A'.-H. Hoffmann et al.
Identification number iterations
of "KFI"
(maximum
ofl number ofl time[ CG-itera-I step[ tions [ size
of perturbation
=
1.0
[cm])
:
I
"true"
I
estimated
I output-
I I
parametervalues
I
parametervalues
I error I
[
equationerror
0
0
I.OOE+03
1 2 3 4 5 6 7 8 9 10 ii
1 3 7 2 7 1 2 7 3 1 1
80000E-03 60000E-04 10000E-04 90000E-03 10000E-05 10000E-03 50000E-03 20000E-05 60000E-04 80000E-03 80000E-03
1.00000E-03 100000E-03 100000E-03 100000E-03 100000E-03 i O0000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03
O.O0000E+O0
6.69672E-02
18
0
1.00E+15
1 2 3 4 5 6 7 8 9 10 11
1 3 7 2 7 1 2 7 3 1 1
80000E-03 60000E-04 10000E-04 90000E-03 10000E-05 10000E-03 50000E-03 20000E-05 60000E-04 80000E-03 80000E-03
1.96451E-03 3.68574E-04 7.45696E-04 3.22246E-03 6.56814E-05 1.17401E-03 2.88081E-03 1.69010E-04 3.72298E-04 1.94215E-03 1.33852E-03
1.87564E+01
1.63545E-02
TABLE 3 Identification number iterations
of "KFI"
(maximum of perturbation
of] number ofl time~ CG-itera-I step] tions ] size
I ~ I
"true" parametervalues
0
0
1.00E+03
1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11.
1.80000E-03 3.60000E-04 7.10000E-04 2.90000E-03 7.10000E-05 1.10000E-03 2.50000E-03 7.20000E-05 3.60000E-04 1.80000E-03 1.80000E-03
20
0
1.00E+15
I. 2. 3. 4. 5. 6 7 8 9 10 11
1 3 7 2 7 1 2 7 3 1 1
80000E-03 60000E-04 10000E-04 90000E-03 10000E-05 10000E-03 SO000E-03 20000E-05 60000E-04 80000E-03 80000E-03
I I I
: 10.0
estimated parametervalues
[cm3)
:
] output~ error ]
equationerror
1.00000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03 100000E-03
O.O0000E+O0
1.82471E-01
6.77973E-04 1.95818E-04 3.23896E-04
1.71934E+02
1.74736E-01
5.40S46E-04 -6.06988E-05 6.94250E-04 -1.65050E-03 7.60959E-04 2.86487E-04 1.29586E-03 -7.62096E-03
Adv. Water Resources, 1991, Vol. 14, No. 5
235
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann e t a [ .
TABLE 4 Identification of "KV" ( maximum of perturbation = number ofl number ofl timeiteraI CG-itera-I steptions I tions I size
I
"true" I parameter- I values I
estimated parametervalues
I output-
[ equationI error I error i upper/lowerl upper/lower
0
0
1.00E+03 1. 2.60000E-09 2. 5.30000E-I0
I.O0000E-09 I.O0000E-09
O.O0000E+O0 O.O0000E+O0
4.71841E-04 4.71841E-04
1
4
1.00E+03 1. 2.60000E-09 2. 5.30000E-I0
2.59430E-09 5.42554E-I0
8.58445E-04 5.52064E-04
1.68795E-06 1.68795E-06
2
7
1.00E+03 1. 2.60000E-09 2. 5.30000E-10
2.59998E-09 5.30287E-I0
4.69683E-06 4.68334E-06
8.40697E-09 8.40696E-09
3
7
1.00E+04 1. 2.60000E-09 2. 5.30000E-10
2.60000E-09 5.30001E-I0
9.76111E-09 9.74435E-09
1.79345E-11 1.79150E-11
4
4
1.00E+OS
1. 2.60000E-09 2. 5.30000E-10
2.60000E-09 5.30000E-10
8.56385E-10 7.64369E-10
7.79037E-13 5.59382E-13
5
0
1.00E+06 I. 2.60000E-09 2. 5.30000E-10
2.60000E-09 5.30000E-10
8.56385E-I0 7.64369E-10
7.79037E-13 5.59382E-13
Identification of "KV"
(maximum of perturbation =
number ofl number ofl timeitera[ CG-itera-I steptions I tions I size
236
I I
0.0 [cm] ):
l
I I
"true" I parameter- [ values I
0.i [cm]) :
estimated I outputI equationparameter- I error I error values i upper/lowerl upper/lower
0
0
1.00E÷03 1. 2.60000E-09 2. 5.30000E-10
I.O0000E-09 1.00000E-09
O.O0000E+O0 O.O0000E+O0
1.78945E-03 1.30531E-03
1
8
1.00E+03 1. 2.60000E-09 2. 5.30000E-I0
2.47581E-09 1.52009E-09
1.96427E+00 1.66552E+00
1.72622E-03 1.22716E-03
2
7
1.00E+03 i. 2.60000E-09 2. 5.30000E-10
2.48093E-09 1.53078E-09
1.96602E+00 1.66708E+00
1.72618E-03 1.22715E-03
3
7
I.OOE+04 1. 2.60000E-09 2. 5.30000E-10
2.48094E-09 1.53101E-09
1.96602E+00 1.66708E+00
1.72618E-03 1.22715E-03
4
4
1.00E+05 1. 2.60000E-09 2. 5.30000E-10
2.48094E-09 1.53101E-09
1.96602E+00 1.66708E+00
1.72618E-03 1.22715E-03
5
0
I.OOE+06 1. 2.60000E-09 2 . 5.30000E-10
2.48094E-09 1.53101E-09
1.96602E+00 1.66708E+00
1.72618E-03 1.22715E-03
Adv. Water Resources, 1991, Vol. 14, No. 5
Adaptive methods for parameter identification in ground water hydrology." K.-H. Hoffmann et al,
TABLE 5 Identification of "KV" (maximum of perturbation = number of[ number of] timeiteraCG-itera-I steptions I tions I size
I I [
"true" parametervalues
I I I
1.0 [cm])
estimated parametervalues
:
] outputI equation] error [ error I upper/lower[ upper/lower
0
0
1.00E+03
1. 2.60000E-09 2. 5.30000E-10
1.00000E-09 1.00000E-09
O.O0000E+O0 O.O0000E+O0
1.72705E-02 1.22704E-02
1
8
1.00E÷03 1. 2.60000E-09 2. 5.30000E-10
1.38415E-09 1.02796E-08
1.96382E+01 1.66552E+01
1.72623E-02 1.22714E-02
2
5
I.OOE+03 1. 2.60000E-09 2. 5.30000E-10
1.38407E-09 1.04961E-08
1.96557E+01 1.66708E+01
1.72622E-02 1.22715E-02
3
7
1.00E+04
1. 2.60000E-09 2. 5.30000E-10
1.38405E-09 1.05010E-08
1.96557E+01 1.66708E+01
1.72622E-02 1.22715E-02
4
4
1.00E+OS
1. 2.60000E-09 2. 5.30000E-10
1.38405E-09 1.05010E-08
1.96557E+01 1.66708E+01
1.72622E-02 1.22715E-02
5
4
1.00E+06
1. 2.60000E-09 2. 5.30000E-10
1.38405E-09 1.05010E-08
1.96557E+01 1.66708E+01
1.72622E-02 1.22715E-02
6
0
1.00E+07
1. 2.60000E-09 2. 5.30000E-10
1.38405E-09 1.05010E-08
1.96557E+01 1.66708E+01
1.72622E-02 1.22715E-02
Identification of "KV"
(maximum of perturbation = I0.0
number ofl number ofl timeitera[ CG-itera-] steptions [ tions I size
I I I
"true" parametervalues
I I I
[cm])
estimated parametervalues
:
I outputI equationI error I error I upper/lowerl upper/lower
0
0
1.00E+03 I. 2.60000E-09 2. 5.30000E-10
1.00000E-09 1.00000E-09
O.O0000E+O0 O.O0000E+O0
1.72664E-01 1.22704E-01
1
13
I.OOE+03 i. 2.60000E-09 2. 5.30000E-10
-I.17109E-08 9.37074E-08
1.93822E+02 1.66518E+02
1.72655E-01 1.22707E-01
2
8
1.00E+03 1. 2.60000E-09 2. 5.30000E-I0
-1.17700E-08 9.58874E-08
1.93993E+02 1.66674E+02
1.72655E-01 1.22708E-01
3
7
1.00E+04 1. 2.60000E-09 2. 5.30000E-I0
-1.17704E-08 9.59373E-08
1.93993E+02 1.66674E+02
1.72655E-01 1.22708E-01
4
9
1.00E+05 i. 2.60000E-09 2. 5.30000E-I0
-1.17704E-08 9.59374E-08
1.93993E+02 1.66674E+02
1.72655E-01 1.22708E-01
5
4
1.00E+06
I. 2.60000E-09 2. 5.30000E-10
-1.17704E-08 9.59374E-08
1.93993E+02 1.66674E+02
1.72655E-01 1.22708E-01
6
0
1.00E+07 1. 2.60000E-09 2. 5.30000E-10
-1.17704E-08 9.59374E-08
1.93993E+02 1.66674E+02
1.72655E-01 1.22708E-01
Adv. Water Resources, 1991, Vol. 14, No. 5
237
Adapttt, e methods for parameter identification in ground water hydrology: K.-H. Hoffmann et al.
TABLE 6 Identification of "KAKST-WATER" number of[ number ofl timeitera] CG-itera-I steptions I tions ] size
( maximum of perturbation =
] I I
[ I I
estimated parametervalues
[cm] ): equationerror
I output] error
I
I.O0000E-09 1.00000E-09 1.00000E-09 1.00000E-09
O.O0000E+O0
1.41416E-03
4.00000E-09 1.00000E-09 5.00000E-10 B.OOOOOE-IO
3.99674E-09 1.00005E-09 5.00567E-10 8.00235E-10
6.69588E-04
1.53712E-06
I.OOE+03 i. 2. 3. 4.
4.00000E-09 1.00000E-09 5.O0000E-IO 8.00000E-10
4.00000E-09 I.O0000E-09 5.00000E-IO 8.00000E-IO
7.38529E-07
1.83602E-10
4
1.00E+04 1. 2. 3. 4.
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-10
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-10
5.37482E-08
7.35190E-13
0
1.00E+05 I. 2. 3. 4.
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-I0
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-10
5.37482E-08
7.35190E-13
1. 4.00000E-09 2. 1.00000E-09 3. 5.00000E-I0 4. 8.00000E-10
0
1.00E+03
8
1.00E+03 1. 2. 3. 4.
11
Identification of "KAKST-WATER" number of] number of[ timeitera ~ I CG-itera-I steptions ] tions J size
238
"true" parameterval~es
0.0
] I
(maximum of perturbation = "true" parametervalues
I I
]
estimated parametervalues
0.1
[cm])
:
[ output] error I
equationerror
0
0
1.00E+03 1. 2. 3. 4.
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-10
1.00000E-09 1.00000E-09 1.00000E-09 1.00000E-09
O.O0000E+O0
1.88894E-03
1
9
I.OOE+03 i. 2. 3. 4.
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-I0
4.03760E-09 9.78280E-10 4.82959E-10 7.95862E-10
1.66541E+00
1.22701E-03
2
9
I.OOE+03 i. 2. 3. 4.
4.00000E-09 I.O0000E-09 5.00000E-10 8.00000E-10
4.04090E-09 9.78202E-I0 4.82370E-10 7.95622E-10
1.66697E+00
1.22701E-03
3
9
1.00E+04
1. 4.00000E-09 2. 1.00000E-09 3. 5.00000E-10 4. B.OOOOOE-IO
4.04090E-09 9.78202E-10 4.82369E-I0 7.95622E-10
1.66697E+00
1.22701E-03
4
0
1.00E+05 1. 4.0bOOOE-09 2. 1.00000E-09 3. 5.00000E-10 4. 8.00000E-10
4.04090E-09 9.78202E-10 4.82369E-10 7.95622E-10
1.66697E+00
1.22701E-03
Adv. Water Resources, 1991, Vol. 14, No. 5
Adaptive methods for parameter identification in ground water hydrology: K.-H. Hoffmann e t al.
TABLE 7 Identification of "KARST-WATER" number ofl number ofl timeiteraI CG-itera-I steptions I tions ] size 0
(maximum of perturbation =
I I I
"true" parametervalues
I ] I
estimated parametervalues
1.0
[cm])
:
I outputI error I
I equationI error I
1.00E+03 1. 4.00000E-09 2. 1.00000E-09 3. 5.00000E-10 4. 8.00000E-I0
1.00000E-09 1.00000E-09 1.00000E-09 1.00000E-09
O.O0000E+O0
1.23779E-02
12
1.00E+03 1. 2. 3. 4.
4.00000E-09 1.00000E-09 5.00000E-10 8.00000E-10
4.40475E-09 7.80676E-10 3.24662E-I0 7.57020E-10
1.66541E+01
1.22701E-02
9
1.00E÷03 1. 2. 3. 4.
4.00000E-09 1.00000E-09 5.00000E-I0 8.00000E-10
4.40837E-09 7.80326E-10 3.23867E-I0 7.56746E-10
1.66697E+01
1.22701E-02
3
9
I.OOE+04
I. 2. 3. 4.
4.00000E-09 I.O0000E-09 5.00000E-IO 8.00000E-I0
4.40837E-09 7.80325E-10 3.23867E-10 7.56746E-I0
1.66697E+01
1.22701E-02
4
0
I.OOE+05 I. 2. 3. 4.
4.00000E-09 I.O0000E-09 5.00000E-IO 8.00000E-IO
4.40837E-09 7.80325E-10 3.23867E-I0 7.56746E-10
1.66697E+01
1.22701E-02
Identification of "KAKST-WATER" number ofl number ofl timeiteraI CG-itera-] steptions ~ tione [ size
(maximum of perturbation
: i0.0 [cm])
:
I [
"true" parameter-
I
[
estimated parameter-
[ outputI error
I equationI error
I
values
I
values
I
I
0
0
1.00E+03
1. 4 . 0 0 0 0 0 E - 0 9 2. 1 . 0 0 0 0 0 E - 0 9 3. 5.00000E-10 4. 8.00000E-10
1.00000E-09 1.00000E-09 1.00000E-09 1.00000E-09
1
12
1.00E+03
1. 4 . 0 0 0 0 0 E - 0 9 2. 1 . 0 0 0 0 0 E - 0 9 3. 5 . 0 0 0 0 0 E - 1 0 4. 8.00000E-10
2
9
1.00E+03
3
9
4
0
O.O0000E+O0
1.22748E-01
7.97770E-09 -1.36264E-09 -1.24130E-09 4.20275E-10
1.66518E+02
1.22702E-01
1. 4 . 0 0 0 0 0 E - 0 9 2. 1 . 0 0 0 0 0 E - 0 9 3. 5 . 0 0 0 0 0 E - 1 0 4. 8.00000E-10
7.98436E-09 -1.36600E-09 -1.24412E-09 4.19760E-10
1.66674E+02
1.22702E-01
1.00E+04
1. 4 . 0 0 0 0 0 E - 0 9 2. 1.00000E-09 3. 5 . 0 0 0 0 0 E - 1 0 4. 8.00000E-10
7.98436E-09 -1.36601E-09 -1.24412E-09 4.19761E-10
1.66674E+02
1.22702E-01
1.00E+05
1. 4 . 0 0 0 0 0 E - 0 9 2. 1 . O 0 0 0 0 E - 0 9 3. 5 . 0 0 0 0 0 E - 1 0 4. 8.00000E-10
7.98436E-09 -1.36601E-09 -1.24412E-09 4.19761E-10
1.66674E+02
1.22702E-01
Adv. Water Resources, 1991. Vol. 14, No. 5
239