Applications of Kalman filtering in the analysis and design of groundwater monitoring networks

in the

analysis and design of groundwater

monitoring networks

F.C. van Geer

TR dissl

in the

in the

analysis and design of groundwater

monitoring networks

Proefschrift

ter verkrijging van de graad van doctor aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus, Prof.dr. J . M . Dirken, in het openbaar te verdedigen ten overstaan van

een commissie aangewezen door het College van Dekanen

op dinsdag 23 juni 1987 te 14.00 uur door

Francois Cornells van Geer

geboren te Naarden Civiel-ingenieur

Prof.dr. P.A. Burrough Prof .dr.ir. J . J . Kalker

ACKNOWLEDGEMENTS

1. INTRODUCTION

1.1 Objective of the study

1.2 Groundwater monitoring networks and data collection 1.3 The quality of a groundwater monitoring network 1.4 Methodology

2. GROUNDWATER FLOW IN TERMS OF SYSTEMS THEORY 2.1 Introduction

2.2 Groundwater flow 2.2.1 Darcy's law

2.2.2 Equation of continuity 2.2.3 Model equations

2.2.4 Finite difference approximation 2.3 The concept of state description

2.3.1 The system

2.3.2 The state of the system

2.3.3 Application of the system definitions to groundwater flow

2.3.4 The state equation 2.4 Measurement equation 2.5 Recapitulation

3. THE OPTIMIZATION OF MONITORING NETWORKS 3.1 Introduction

3.2 Optimal design of groundwater monitoring networks 3.3 Groundwater monitoring network analysis

3.4 Estimation principles

3.5 Recent developments in hydrological monitoring network design

3.5.1 Time series analysis and regression methods 3.5.2 Kriging

3.5.3 Kalman filtering 3.6 Choice of Kalman filtering

4.1 Introduction 55 4.2 The calculation scheme of the Kalman filter algorithm 55

4.2.1 Estimation properties 55 4.2.2 Starting points 56 4.2.3 Calculation scheme of the Kalman filter algorithm 57

4.3 Application to groundwater flow 63

4.3.1 The input vector 64 4.3.2 Initial conditions 64 4.3.3 Measurements and related quantities 65

4.3.4 The state equation and related quantities 66

4.4 The use of the filter algorithm 71

4.4.1 Interpolation 72 4.4.2 Simulation 74 4.4.3 Forecasting 75

CASE STUDY: APPLICATION TO THE RUURLO AREA 77

5.1 Introduction 77 5.2 A brief description of the geohydrology 78

5.3 Modelling 80 5.4 Calibration 85 5.5 Calculation results 85

5.5.1 Monitoring network density 85 5.5.2 Trade-off between network density and measurement

frequency 91 5.6 Discussion 93

DISCUSSION 97 6.1 Introduction 97 6.2 The influence of modelling on the estimation errors 97

6.2.1 The measurement equipment 98 6.2.2 The network density and the network frequency 98

6.2.3 The choice of the model 98 6.3 Deterministic and stochastic models 99

6.3.1 Deterministic parameters 100 6.3.2 Stochastic parameters 101 6.3.3 Comparison between deterministic and stochastic

1. In de opleiding van de Nederlandse hydroloog krijgt de numerieke wiskunde te weinig aandacht om geohydrologische modelstudies verantwoord uit te voeren.

2. Bij de modellering van grondwaterstroming op basis van hetzij een zuivere deterministische beschrijving danwei een zuivere stochastische be­ schrijving, wordt niet ten volle gebruik gemaakt van de in het waarnemings­ materiaal aanwezige informatie.

3. De zwakste schakel bij het ontwerpen van een grondwaterstandsmeetnet is het exact formuleren van de doelstelling ervan.

4. Het optimale grondwaterstandsmeetnet bestaat niet zelfs indien de doel­ stelling exact is geformuleerd en kan dus niet ontworpen worden.

5. Om een betere evaluatie van de gevolgen van een grondwaterwinning mogelijk te maken, zou behalve een meetnet ook debiets fluctuaties bij de vergun­ ningverlening verplicht gesteld moeten worden.

6. Bij verlagingsberekeningen is tijdreeksanalyse superieur aan regressie analyse. (Box, G.E.P. and Jenkins, G.M., 1960, Time series analysis: forecasting and control. Holden-Day, San Francisco.)

7. De modelcoëfficiënten van een tijdreeksmodel voor een grondwaterstandsreeks zijn af te leiden uit de geohydrologische parameters van het desbetreffende gebied.

8. Het verdient aanbeveling een nauwere samenwerking tot stand te brengen tussen instellingen die zich bezighouden met stochastisch modelleren van hydrologische processen. (Grondwatermeetnetten en Databestanden, Verslag van een workshop 17-18 april 1986. Rapportnummer DGV-TNO: PN 86-11).

door een indicatie van de betrouwbaarheid aan te geven. (Marsman en De Gruiter, 1984. Dutch soil survey goes into quality control. In: Soil information systems technology (ed. P.A. Burrough and S.W. Bie).

PUDOC, Wageningen.)

10. Veel toegepast wiskundigen onderschatten de problemen bij het toepassen der wiskunde.

F.C. van Geer. Applications of Kalman filtering in the analysis and design of groundwater monitoring networks.

6.5 Interpretation 106 6.5.1 Spatial interpretation 108

6.5.2 Interpretation in respect to time 109 6.6 Groundwater monitoring network design 110

6.6.1 The criterion 110 6.6.2 The measurements 112 6.7 Extension of the application of Kalman filtering 113

6.7.1 Control problems 114 6.7.2 Non-linear problems 115 6.7.3 Recommendations 118 REFERENCES SUMMARY SAMENVATTING APPENDICES

I DETERMINATION OF THE MEASUREMENT MATRIX

II DERIVATION OF THE KALMAN FILTER ALGORITHM

III INDEPENDENCE OF THE INNOVATIONS

IV THE CISKA COMPUTER PROGRAM PACKAGE

This thesis results from the CT-7904 Research Project of the Delft University of Technology. The research work was carried out at the Discipline Group of Water Management of the Department of Civil Engeneering, in cooperation with the TNO-DGV Institute of Applied Geo­ science.

I am very grateful to the project leader Ir. P. van der Kloet for introducing me to and guiding me through the labyrinth of stochastic modelling.

I also wish to thank the staff of the Monitoring Department of the TNO-DGV Institute of Applied Geoscience for geohydrological advice and their support in the realization of the computer programs. In particular, I wish to thank Ms. M.C.A. Cornel-Reniers and Ms. A.A.C. Gieskens for typing the text and the numerous formulae and Ms. S.M. Ploos van Amstel for preparing all the figures.

in the

Chapter 1 INTRODUCTION

1.1 Objective of the study

Groundwater plays an important role in the basic needs of human society, in industrialized as well as developing countries. It is both a raw material for water supply and agriculture and a load on subsoil constructions and dams. Often insight is needed into the natural fluctuations of the groundwater level, changes in piezometric levels (defined in Chapter 2) and groundwater velocities due to human activities and also into the prediction of changes due to activities proposed. This insight is gained by measurement of the actual behaviour of groundwater and by theoretical knowledge about the physical laws describing the process of groundwater movement. Measurements of piezometric levels are taken in observation wells. If a number of observation wells are installed for the same objective, those we]Is are called a groundwater monitoring network.

At first, observation wells were placed at more or less arbitrary locations. Until recently, little research had been done on the quality of groundwater monitoring networks. A few years ago, however the DGV-TNO Institute of Applied Geoscience started the optimization of the primary groundwater monitoring networks in the Netherlands. The optimization procedure consists of two steps, each performed separately:

1. Optimization of the network density.

2. Optimization of the measurement frequency, where only observation wells are used, at locations resulting from Step 1.

Both steps are carried out by purely statistical methods; the knowledge about the groundwater system is used indirectly.

The objective of the present study is:

To develop a method for analysing the quality of existing groundwater monitoring networks, and to give guidelines for the optimization of those networks, in such a way that physical knowledge about the groundwater system and the relation between network density and measurement frequency are included in the method.

The method of data collection and the types of groundwater monitoring networks are discussed in Section 1.2. An explanation is then given of what is meant by the quality of a monitoring network in Section 1.3. Finally in Section 1.4 the methodology used in the study is briefly introduced.

1.2 Groundwater monitoring networks and data collection

Most groundwater measurements are taken to obtain information to describe the quantitative aspects of groundwater flow. These quantitative aspects include elevation of groundwater tables, piezometric levels, groundwater velocities, groundwater discharges and volumes of available groundwater. All quantitative aspects are detectable from measurements of the piezometric level, if some knowledge about the geohydrological conditions is available (see Chapter 2 ) . Therefore a network which is designed to monitor the quantitative aspects of groundwater flow can be restricted to a monitoring network for piezometric levels.

In the last ten years there has been increasing interest in the quality of groundwater. To obtain information about chemical groundwater quality, a sample can be taken from the water in an observation well and analysed. The analysis of water samples is time consuming and expensive. Moreover, most groundwater monitoring networks are designed to gain insight into the quantitative behaviour of groundwater movement and, generally, such a network will not give sufficiently representative information about the extent and changes in groundwater quality processes. Throughout this thesis only those monitoring networks are discussed, which are designed to observe groundwater flow quantitatively.

Because the process of groundwater movement takes place beneath the land surface, it cannot be observed directly. Therefore to obtain information about the process of groundwater movement, observation wells are made. Such an observation well consists of a drilled hole in which a measurement filter is placed. The water level in the observation well is assumed to be equal to the piezometric level in the aquifer at the location of the measurement filter. The water level in the observation well can be measured very easily.

In the Netherlands, a groundwater monitoring network which covers the whole country is managed by DGV-TNO Institute of Applied Geoscience. This monitoring network consists of many thousands of observation wells. The majority of these wells are monitored either with a frequency of 24 times/year or with a frequency of 4 times/year, resulting in an enormous amount of data. The data storage and data management is an integral part of the monitoring network.

Groundwater velocities and piezometric levels at times and places other than the monitoring times and places are estimated in some way or another. Estimation requires a model (not neccessarily mathematical) in order to relate the observed levels to the levels to be estimated.

1.3 The quality of a groundwater monitoring network

In principle, a groundwater monitoring network is designed to serve some predefined goal. This goal can be the monitoring of the annual fluctuations of the groundwater level in an agricultural area or the changes in the piezometric level in the area surrounding a pumping station. Groundwater movement is a relatively slow process in comparison to surface water flow, measurement series of several years or even several decades are needed to provide sufficient information. In most cases, however, a problem has to be solved in a short period of time and hence there is no time to collect long measurement series. To overcome this problem multi-purpose groundwater monitoring networks have been designed. At the design stage of such networks often the applications of the measurements are not known.

Groundwater monitoring networks can be classified by means of the use which is or will be made of the measurements. For example, distinction can be made between regional and local networks or between permanent and temporary networks. Multi-purpose monitoring networks are permanent and nearly always regional networks, while monitoring networks designed for research objectives are generally local and temporary. Measurements from different networks may be used for different purposes, but from the point of view of monitoring network management, they can be dealt with in the same way. Therefore, although the present study is focussed on the groundwater monitoring

network of the DGV-TNO Institute of Applied Geoscience, the results of the study are applicable to other groundwater monitoring networks. As stated above, groundwater measurements are used to obtain information- about the behaviour of the quantitative aspects of groundwater flow. The measurements are used to estimate the relevant variables (e.g. plezometerlc levels). It will be obvious that, in general, the accuracy of the estimates is related to the amount of information available and there is a positive correlation between the accuracy of the estimates and the amount of information available. Of course the way in which the estimates are made also influences the accuracy of the estimates (see Chapter 3 ) .

The quality of a groundwater monitoring network is generally recognized by hydrologists to be:

The accuracy with which the geohydrological variables can be estimated from measurements of the monitoring network.

The standard deviation of the estimation error has been chosen as a criterion for the quality of a groundwater monitoring network (WMO 1981). The estimation error and the estimation error standard deviation are illustrated in Figure 1.1. Close to the measurements the estimation standard deviation is small, in contrast in between the measurements it is large.

*

/ /

/ /

/ /

N s

^ - 1 ^ " "

:

- - - J

r~

\ \\

—.•*-"'" ^

i-W;

' *

\ \ > * - - true piezometric level \ \ / \ \ - estimated piezometric level

: V ** \ v- estimation error

v— -N estimation error standard deviation

■ x (distance)

1.4 Methodology

To achieve the objective of the present study a method is needed to calculate estimates of piezometric levels and the error standard deviations of those estimates. Up to the present, few attempts have been made in the optimization of geohydrologlcal monitoring networks to develop a method that takes into account the spatial and temporal behaviour of groundwater and uses available physical knowledge about the groundwater system at the same time.

Analogous problems occur in fields other than hydrology. The problems can be generalized into abstract mathematical formulations and general solutions can be obtained. This is the field of systems theory. In systems theory, an advanced method is available which satisfies all demands of the objective mentioned above. This method is called Kalman filtering. However, a number of difficulties have to be overcome before Kalman filtering can be applied to problems in groundwater monitoring network optimization. The main part of this thesis deals with the adaptation of Kalman filtering in such a way that it is applicable to groundwater flow.

The most important conclusion of this study is that Kalman filtering is a very powerful method for the analysis of groundwater monitoring networks. Therefore it should become a basic tool in groundwater monitoring network optimization.

Chapter 2

GROUNDWATER FLOW IN TERMS OF SYSTEMS THEORY

2.1 Introduction

The aim of this chapter is to present to the reader the concept of systems theory applied to groundwater flow. It will be shown that all elements of the deterministic description of groundwater flow can be expressed in terms of systems theory. These terms will be used in the remainder of this thesis. As an example of the deterministic description of groundwater flow a finite difference model is worked out. This model has been used in the case study discussed in Chapter 5.

The objective of a groundwater monitoring network is to investigate the behaviour of the groundwater in time and space. The goal of the present study is to develop a method for analysing the quality of a groundwater monitoring network (see Chapter 1 ) . As discussed in Chapter 3, the standard deviations of the estimates of the piezometric level are used as a measure of the quality of the network. A Kalman filter algorithm has been chosen for calculating the estimation error standard deviations. The choice of Kalman filtering, which is a well known method in systems theory, is discussed in Chapter 3. The Kalman filter algorithm is applicable for any system which can be described by:

- a linear state equation. This (vectorial) difference equation describes the behaviour of the system;

- a linear measurement equation, which describes the relationship between measurements and the corresponding variables of the system. In this chapter, these two equations are derived for groundwater flow. The state equation can be based on any linear mathematical model. The model can be either deterministic or stochastic. The advantages and disadvantages of both types of models are discussed in Chapter 6. Section 2.2 presents the groundwater flow model chosen to describe the geohydrological conditions. It is shown that the model can be written as a vectorial linear difference equation by means of a finite difference approximation.

Then, in Section 2.3, some basic elements of systems theory are defined. The geohydrological interpretation of these elements is given

and the state equation based on the selected model is derived. In Section 2.4 the measurement equation is derived for measurements of the piezometric level. A recapitulation of the most important results is given in Section 2.5.

2.2 Groundwater flow

The state equation discussed in Section 2.3 is one of the basic equations of the filter algorithm derived in Chapter 4. For groundwater flow, this state equation can be based on any linear mathematical model which is discretlzed in time and space. This means

that the common groundwater flow models (finite elements, boundary elements and finite differences) can be used. A finite difference model is used in the computer program developed. Because the goal of this study is the optimization of regional groundwater monitoring networks (see Chapter 1 ) , only large scale models have been considered. The concept of aquifers and semipervious layers has therefore been adopted (see Bear 1979). The following assumptions have been made:

- The vertical gradient of the piezometric level of the groundwater inside the aquifers is negligible.

- The groundwater itself is homogeneous. This means that the density of the groundwater is independent of time and space.

- Each aquifer is isotropic.

2.2.1 Darcy^s_law

The basic equation with which to describe groundwater flow is Darcy's law, which states that the specific discharge in the s-direction is proportional to the derivative of the piezometric level in the s-direction, that is

k Ü (2.1)

is the specific discharge in the s-direction [L.T ] is the permeability [L.T~ ]

is the piezometric level [L], which is defined by: qs =

where: q ^s k

= p/pg + z (2.2)

where: p is the pressure with respect to a reference —1 -2 pressure (atmospheric pressure) [M.L ,T 1 p is the density of the groundwater [M.L ] g is the acceleration of gravity [L.T ]

z is the elevation with respect to a reference level [L]

2.2.2 £9.uation_of continuity

The second equation of importance when describing groundwater flow is the equation of continuity. This equation can be derived by considering the water balance of a column (dx.dy) of an aquifer. Two types of aquifers should be distinguished.

1. Aquifers with semi-confined groundwater, which are bounded by semi-pervious or impervious layers at the upper and the lower side (see Figure 2.1). An aquifer with semi-confined groundwater is always saturated with groundwater.

|1«*d<1x

Figure 2.1: Column dx.dy of an aquifer with semi-confined groundwater.

2. Aquifers with phreatic groundwater, which are not bounded by any semi-pervious or impervious layer at the upper side (see Figure 2.2). Phreatic groundwater has a free groundwater table.

Figure 2.2: Column dx,dy of an aquifer with phreatic groundwater.

The equation of continuity has been derived for semi-confined groundwater first. Since there is no loss of groundwater the difference between the inflow and the outflow of the column (dx.dy) during an infinitesimal lapse of time dt equals the amount of groundwater stored in the column. It can be shown (Bear 1979) that the storage during dt, with a good approximation, is proportional to the rise of the pressure during dt. Because the groundwater is homogeneous the storage is also proportional to the rise of the plezometric level

(see Eq. 2.2.) so it yields:

D D

- * dq dz.dydt - J dq dz.dxdt + (q.-q-) dxdydt = Sd<|>dxdy

where: S is the storage coefficient of (semi)confined groundwater [-]

D is the thickness of the aquifer [L]

q is the inflow through the upper side [L.T ]

1 i q is the outflow through the lower side [L.T ]

2 i q is the specific discharge in x-direction [L.T ]

x i q is the specific discharge in y-direction [L.T ]

D 3qv - ' — dz o 3x D 3q, _ ; __; o 3y

Division of Equation (2.3) by (dx.dy.dt), with dx, dy and dt tending to zero, gives the equation of continuity for semi-confined groundwater:

n an

dz - q

+ q

= + S U (2.4)

The derivation of the equation of continuity for phreatic groundwater is slightly different from that for semi-confined groundwater. In the case of phreatic groundwater the process of groundwater storage is dominated by the change of the thickness of the saturated zone (h). The storage is proportional to the rise of the groundwater table. Because of the absence of a semi-pervious layer above the aquifer the term q. of Equation (2.3) is absent. Instead, there will be recharge (positive or negative) of the groundwater due to precipitation and evapotranspiration, so:

- J dq dz.dydt - ' dq dz.dxdt + (N-q„) dxdydt = ydhdxdy (2.5)

where: N is the recharge as a depth per unit of time [L.T ] u is the storage coefficient of phreatic groundwater [-] h is the thickness of the saturated zone [L]

Division of Equation (2.5) by (dx.dy.dt), with dx, dy and dt tending to zero, gives the equation of continuity for phreatic groundwater:

h 3q h 3q

/T - ^ d z - ' r1 dz - q, + N = + ii ^ (2.6)

o 3x o 3y H2 M 3t

2.2.3 Model equations

The model equations for semi-confined and phreatic groundwater are ob­ tained by combining Darcy's law and the equation of continuity.

The model equation for semi-confined groundwater has been derived first. By substitution of Equation (2.1) in the first two terms of Equation (2.4) the integrals in Equation (2.4) can be expressed in piezometric levels. These integrals are then approximated by:

ofe'-^^-^H*

' I- {- k !£} dz

*

-

Df-o 3y 3y 3y 3y t K 3y ' (2.8)

where the quantities -r— {k -r-"-} and -r— {k -r^} are assumed to be average values over the interval (0,D) in the z-direction.

The inflow and outflow through semi-pervious layers at the upper and the lower boundaries of the aquifer can be expressed in the piezometric levels in the aquifers. The vertical inflow q1 through a

semipervious layer with thickness d1 and coefficient of permeability

k. (see Figure 2.3) can be formulated according to Darcy's law as:

(2.9)

where: $ i s the piezometric level above the semi-pervious layer [L]

c ^ d . / k . i s the hydraulic resistance of the semi-pervious layer [T] aquifer semi-pervious layer aquifer * i

?m7/zzm

»z

Z P

By analogy the vertical outflow q. may be formulated by:

where: <j>0 is the piezometric level below the semi-pervious layer [L]

c9=d./k„ is the hydraulic resistance of the semi-pervious

layer [T]

Substitution of Equations (2.7), (2.8), (2.9) and (2.10) in Equation (2.4) gives the model equation for semi-confined groundwater:

D f- {kf*} + D |- {k|*} - {-L + -L} $ + -L + -2 = S |i (2.11)

3x 3x 3y 3y c, c' y c1 c„ 3t

In the case of a phreatic aquifer the upper boundaries of the integrals in Equation (2.6) are functions of time. Fortunately in a broad class of problems, it is sufficient to approximate the height of the upper boundaries with the time average of the thickness of the saturated zone H. This will be a good approximation if the fluctuations of the groundwater level in time are small in respect to H, which still may be a function of x and y (Bear 1979). Equation (2.1) is substituted in the first two terms of Equation (2.6). Similarly to the integrals in Equation (2.4) the integrals in Equation (2.6) are approximated by:

Vi H n 1 . ^ d z ^ _ H |- {k |4} (2.12) o 3x 3x 3xJ and h 3qv a a* ; _Zclz -v, _ H f- {k |*} (2.13) o 3y 3y 3yJ

The errors due to the use of the average thickness H are respectively:

f- {k eu |4} and |- {k eu |*}

3x h 3x 3y h 3y where: e, = H-h

h

Because the vertical gradient of the piezometric level is neglected (Assumption 1, Section 2.2) the time derivative of the thickness of the saturated zone is equal to the time derivative of the piezometric level. Using this equality and substituting Equations (2.12), (2.13)

and (2.10) in Equation (2.6) gives the model equation for phreatic groundwater:

& <■&>

+

"I? <4f> - è ?

+

s f

+

» - >&

(2

-

14:

If there are a number of aquifers and semi-pervious layers above each other, the behaviour of the groundwater inside the aquifers can be described by a model equation for each aquifer. The model equations for the aquifers are linked by the vertical flow through the semi-pervious layers. The set of model equations must be solved simultaneously.

In the derivation of the model equations a number of assumptions have been made. A summary is given below of the conditions that have to be

fulfilled in order for Equations (2.11) and (2.14) to be valid.

1. The structure of the subsoil should be such that it can be described schematically as a number of horizontal layers

(aquifers, semi-pervious layers and impervious layers) above each other.

2. The vertical gradient of the plezometric level in the aquifers is negligible, which implies that the vertical component of the groundwater flow is small. In large scale problems where the horizontal dimensions are much larger than the vertical dimension, this condition will nearly always be fulfilled. However, in small scale problems, such as groundwater flow through a dike, the model equations may not give a good approximation.

3. The fluctuations in time of the groundwater level in a phreatic aquifer must be small with respect to the time average of the thickness of the saturated zone in that aquifer, at least during the time interval of interest.

4. Each aquifer is assumed to be isotropic. This assumption is made so that the derivation of the model equations will be clear as possible. Model equations for non-isotropic aquifers can be derived in a way similar to that given in this chapter, but, in

large scale problems, there will seldom be sufficient data to fill in the actual parameter values of the non-isotropic aquifer.

2.2.4 Finite difference approximation

In general, the model Equations (2.11) and (2.14) are difficult to solve analytically. For this reason, the model equations are usually solved with the help of numerical methods, where the continuous piezometric surface is calculated at a finite number of points in space and time. The model equations have therefore, to be discretized in space and time. The most common methods of solving the equations numerically are methods based on approximations with finite differences, finite elements or boundary elements. In this study the finite difference approach has been chosen.

Both model Equations (2.11) and (2.14) are differential equations of the diffusion type, which can be approximated by difference equations. Because all parameters (k, D, H, S, y, c. and c?) are positive quantities, the initial conditions will damp out in a finite period of time (Mitchel 1971).

First the space derivatives are approximated by differences and then the time derivative are dealt with.

a) Space derivates

The area under consideration space derivatives is subdivided by a number of gridlines (see Figure 2.4).

The lines i=0, i=n and 1=0, j=n are the boundaries of the area. The

x J J y

gridlines need not to be equidistant. Together the lines form a calculation grid. The space derivatives at all gridpoints (i,j)

(defined in Figure 2.4) can be approximated by the differences, Equations (2.15) and (2.16).

0,n» OJ+1 O.j O.j-1 Ayj A> H A xi - i i-1J Ax-, U*1 i.j U-1 '♦1,j

i-1,0 i,0 i+1,0 nx, 0

F i g u r e 2 . 4 : S u b d i v i s i o n of a n a r e a b y l i n e s in x a n d y d i r e c t i o n .

- i ( ki , j+ ^ - i . j ) (*i,j - * i - i , j) / i xi - i}' iu**±-i+ Axi> > ( 1=1,..., n -1)

(J-0, ny) (2.15)

( 1=0,..., nx)

where ij> . is the plezometric level at grid point (i,j) [Lj Ax is the distance between the lines i and (i+1) [L] Ay is the distance between the lines j and (j+1) [L]

k± is the permeability at grid point (i,j) [L.T-1]

If the gridlines are equidistant with Ax. = Ax and Ay., = Ay. and the aquifer is homogeneous with k = k, the differences, Equations

(2.15) and (2.16) become

f

{&} I . k *i+i.j ' ]*±4

*i-i.3 '

.

3x 3xJ 'i,j ( A x ) *

and

*- { k ^ } I * k "

'J

"

*

'J

^ ' i -

(

16a)

3y

* V 'i.J ( A ? P

'

At grid point (i,j) the model equation for semi-confined groundwater can be approximated by substituting Equations (2.15) and (2.16) in Equation (2.11) giving:

3t 'i.J = al*i,j-l+ a2'',i-l,j+ a3*i,j+ 04*i+l,j+ a5*i,j+l+ Bl*l+ 62*2 ( 2'1 7 )

where: ai = { ( ki , j + ki , j - i) - D } ' ^ y j - i ^ y j +fiyj_1)-s} a2 = { ( ki , j + ki - l , j) - D } / ^xi _ i ^xi + A x ^ p . S } a4 = { ( ki + l , j + ki , j) - D } ' { A Xi ( A Xi +A xi_1) . S } a5 = { ( ki , j + l + ki , j) , D } ' { A yj ( A yj + AyJ_1).S} '1 S .C l 5

= _ L _

The quantities D, S, Cj and c2 are taken at grid point (i,j).

A similar equation is obtained for phreatic groundwater by substituting:

D = H, S = u, ^ = N and Bx = 1/v (2.18)

Groundwater abstraction

Assume that a volume Q is abstracted (or injected) per unit of time i»J

at grid point (i,j). Such an abstraction can be taken into account by adding an input term to the righthand side of Equation (2.17). This can be done by considering the area AR. (see Figure 2.5) around grid

i»J point (i,j).

For confined groundwater the contribution to the rate of change of the piezometrlc level at grid point (i,j), due to the abstraction, can be expressed by:

%L

r

y/z

*\y /

30/

/ / /

>

w.j

The contribution to the rate of change of the piezometric level in grid point (i,j) for phreatic groundwater can be derived from Equation (2.19) by substituting the storage coefficient of phreatic groundwater u for the storage coefficient of confined groundwater S. The model equation for grid point (i,j) in a confined aquifer with an abstraction of groundwater is:

3^ *1.J

at 'i,j" otl*i,j-l+a2*l-l,j+ot3*i,j+ot4(f>i+l,j+a5*i,d+l+6l<*,l+624,2+ S.AR. , i»J

(2.20) With (2.18) a similar model equation can be obtained for phreatic groundwater.

Equation (2.20) can be written in vector notation:

3£

8 t U , J

[°l

2

3

4

5]

+ [8^

To calculate the piezometric level in all grid points, boundary conditions are required.

There are two types of boundary conditions: - The piezometric level at the boundary is given. - The groundwater flux across the boundary is given.

It should be noted that a combination of the two types of boundary conditions is also possible. In the computer program (see Chapter 4) only given piezometric levels are used as boundary conditions.

For ease of manipulation it is preferable to have all known variables in one vector. Therefore the vector Equation (2.21) is modified for the grid points on the lines i=l, i=n -1 and j=l,

j=n - 1 . y

For example, for grid point (i,l) the piezometric level at boundary point (1,0) is a given boundary condition. Therefore Equation (2.21) for grid point (1,1) is modified:

3*1 3t'i,l

K

3V5

•

V

•J

h

For the grid points on the lines i=n - 1 , j = l and j=n -1 analogous modifications are made.

The model equation for each grid point can be written in vector notation, such as Equations (2.21) or (2.22). To solve the set of equations simultaneously, all vectorial model equations are combined to one vectorial equation. Therefore the vectors <f> and u are defined. The vector $ comprises the piezometric levels at all grid points except the boundary points.

(2.23) '1,1 >2.1 nx-1,1 "1,2 nx-l,ny-l

The second vector is the input vector. The physical processes which are represented by the terms <(> , <(> , N, Q and by the boundary

1 ^ 1, j

conditions are quite different from each other. Nevertheless from a mathematical point of view they may be dealt with in the same way. Therefore the input vector u is defined as:

(2.24)

where: u. is a vector whose components are the $ , <)>„ and N values at each grid point.

u„ is a vector whose components are the abstractions of groundwater at grid points where groundwater Is abstracted.

u, is a vector whose components are the piezometric levels at the boundary points.

The lefthand sides of all model equations form the time derivative of the vector <(>. The right hand side of each model equation consists of two terms. The first term contains a vector with the dependent variables (piezometric levels) and the second term contains a vector with the controlling variables (input, abstractions, boundary

conditions). If there are several aquifers above each other, all dependent variables are combined into the vector if> and all controlling variables are combined into the vector u. The set of model equations can be written as one vectorial differential equation using the vectors

<t> a n d u .

3t A c c <(> + B u (2.25)

where: A and B are matrices dependent on the coefficients a. through

a5 and 6 ^ 62> 1/S.A^ and 1/u.AR^ .

Equation (2.25) is the model equation which is discretized in space.

b) Time derivative

The next step is the discretization in time. The approximations, Equations (2.26), (2.27) and (2.28) are therefore made during the timestep from (k-l)At to kAt.

3t <)> <\/ 0 'k ^ At 5* k u "^ u. k where

V

= *

In the following the symbol k is used to indicate the discrete time unless it is clear from the context that the coefficient of permeability is meant. Substitution of Equations (2.26), (2.27) and

(2.28) in (2.25) and rewriting the result yields:

K =

k (2-29)

with A = [i - 0.5 At.A ]- 1 . [i + 0.5 . AtA ]

B = [i - 0.5 At.A l"1 . B

L cJ c

I is the unit matrix.

The vectorial model equation derived is known as a difference equation based on the Crank-Nicolson scheme. Accordingly to Bear (1979) this scheme is convergent and unconditionally stable.

2.3 The concept of state description

2.3.1 The system

There are several definitions in use for the basic elements of systems theory. No attempt is made here to give extensive definitions of those elements. They are defined simply in a way suitable for this study. An extensive discussion about this subject can be found in Kwakernaak & Sivan (1972) and Ogata (1969).

The starting point of an approach using systems theory is reality. In principle reality is unbounded. It is assumed that in reality physical processes occur. Groundwater flow is such a process. In each particular case one is not interested in complete reality or in all physical processes. Depending to the problem which has to be investigated and the relevant processes, part of reality, can be separated from its environment by boundaries. The separated part is called a system. Within the system there are variables, called the system variables. There may be relationship between the system variables. Such relationship are dependent on the properties of the system, characterized by parameters. The environment influences the system by means of input variables, which are independent of the system. In reverse, the environment is influenced by the system by means of output variables.

It is clear that all output variables are system variables. A system with boundaries, variables and parameters is sketched in Figure 2.6.

environment

input variables

Figure 2.6: A system

In this study the following definitions are used. - System

a part of the reality separated from its environment by boundaries.

- System variable

a variable inside the system.

- Parameter

a quantity which characterizes a property of the system.

- input variable

a variable outside the system, which influences the system. SYSTEM system variabtes parameters output variables system boundary

- output variable

a system variable which influences the environment.

In each particular case the system, the system boundaries, the variables and the parameters have to be chosen according to the phenomenon to be studied.

2.3.2 The state of the system

To study the behaviour of the system a representation of the process in the system can be made. Such a representation is called a model. There are several types of models, but in this study only mathematical models are under consideration. A mathematical model of a system is obtained by a schematical and simplified description of the relationship between the system variables, the input variables and the parameters. The relevant relationship are written as an equation or a set of equations. By solving the equations - analytically or numerically - the behaviour of the variables in the model is known and this behaviour can be interpreted to gain insight into the behaviour of the system variables.

A system is referred to as causal if the system variables react only to past and present inputs. Then the system variables, at time t , are completely determined by all input variables at times t (-co < t < t ) . The influence of the input (- <» < t < t ) can be stored in a quantity

o

which is called the state of the system which, in this thesis, is defined as:

The state of the system is a set of values of system variables such that all information of the past, of the system, relevant for the future behaviour, is embedded in this set of values.

According to this definition the future behaviour of the system variables is completely determined by the present state of the system and the future input. As a consequence of the definition of the state, the system variables, which are not part of the state at time t , are determined by the state at that time.

2.3.3 A££lication_of_the_sy_stem_definitions_to groundwater_flow

- System

In problems concerning flow of homogeneous groundwater in the saturated zone a system always contains a separated volume of that zone. The choice where the boundaries of the system are located is dependent on the geographical extent of the area under consideration, the depth at which information is required and on the problems to be solved.

- System variables

There are three types of system variables: . piezometrlc levels

. groundwater velocities . volumes of groundwater.

- Parameters

The parameters are the geohydrological properties of the subsoil.

- Input variables

The input variables are all influences from the environment on the system. The most important input variables are recharge due to precipitation and groundwater abstractions.

- Output variables

The output variables are groundwater levels and discharges.

- State of the system

To define the state of a groundwater flow system it is helpful to know that the velocities and volumes are completely determined by the piezometrlc surface. The future behaviour of the system can be described completely by the present piezometrlc surface and the future input. Therefore the state of the system at time t may be defined by the piezometrlc surface at that time. There are other ways of defining the state of a groundwater system. In fact, there

is an infinite number of possibilities because any linear function of the piezometrlc surface will satisfy the definition of the state of the system. Therefore the state of the system could also have

been defined by using groundwater velocities and a reference level. In the present study, which aims at the optimization of quantitative groundwater monitoring networks, it is most useful to define the state of the system as the piezometric surface.

2.3.4 The state equation

Because the most common groundwater flow models are linear and are modelled discretely in time and space, only the linear discrete form of the state equation is discussed here. Since the state of the system is defined as the piezometric surface, the elements of the vector <j>, form the discretized state of the system. In the remainder of this thesis the vector <j>, is referred to as the state vector.

From systems theory it is known that the evolution of the state of any linear system in discrete time can be modelled by a difference equation of the form of Equation (2.29), see Zadeh & Desoer (1963). The flow diagram of the model is given in Figure 2.7, where <j> is the initial value of the state vector.

to

♦k-1

uk

model ♦k

Figure 2.7: Flow diagram of the model

To build a model, such as Equation (2.29), assumptions and simplifications must be made and as a consequence the piezometric level calculated with this equation will differ from the real piezometric level. The system noise is defined as the difference between the solution of the model and the real state. The state equation, (2.30) can be obtained by adding a vector representing the system noise to the model Equation (2.29).

*k = A *k-l + B Uk+ Wk ( 2-3 0 )

Because the term w, represents all unknown influences which are not modelled in the other terms of Equation (2.30), the vector w, may be considered as a realisation of a stochastic process. A brief introduction to stochastic processes is given below. A detailed discussion about stochastic processes is given by several authors, for example by Papoulis (1965).

Stochastic variable

The basic idea behind stochastic variables is a probability experiment, that can be seen as an action resulting in an a priori unknown outcome. A stochastic variable is defined as a function of the outcome of a probability experiment. The function argument of a stochastic variable is nearly always dropped in the notation. The actual value of a stochastic variable is the result of one particular outcome. This actual value is called the realization of the stochastic variable.

Probability density function

Each single outcome of a probability experiment is a priori unknown, but if the experiment is repeated a number of times, the realization will show some distribution, a probability distribution. The probability distribution can quantifified by the probability density function.

Let the possible realizations of a stochastic variable x be the value a (-•» < a <+<»). The probability that a realization of the stochastic variable x is in the range (-» < a) is a figure between zero and one denoted by:

Px (x < a) (2.31)

The probability density function of x is defined by:

P (x < a + da) - P (x < a) (2.32)

px(a ) . .* _ *

Where: p (a) is the probability density function. da is an infinitesimal range.

The probability density function is scaled such that:

ƒ Px (a) da = 1

(2.33)

EXAMPLE

For a Gaussian variable x the probability density function is given by:

P. (a)

_1Ü

~P

t"

₁

^]

where : u and 0 are the parameters of the probability density function.

The symbol y used here is not related to the storage coefficient of phreatic groundwater.

The graph of the Gaussian density function is given in Figure 2.8.

, P > )

Figure 2.8: Gaussian density function.

Joint probability density function

Often, more than one stochastic variable is of importance. Let the possible realizations of the stochastic variables x.. and x„ be respectively a, and a-. The probability that both the realization of x. is in the infinitesimal range (a. < x < a, + da.) and the realization of x. is in the infinitesimal range (a2 < x < a, + d a2) is

given by:

where: p (a,. a~) is the joint probability density function.

X , t 3C* 1 ^

Expression (2.35) can be generalized for n stochastic variables

x., x„, ...,x with realization a., a„...,a . Analogous to the probability density function of a single stochastic variable, the joint probability density function is scaled such that:

+ 0 0 +0O + 0 0

f f •■• f Px , , x9, . . . , x < V V - - . . °n) d«l dVd Qn - l ( 2-3 6 )

—oo —oo —oo 1 2 n

If the stochastic variables x.,x.,...,x are the components of a stochastic vector x and realizations a.,a.,...,a are the components of the vector a, the joint probability function for x^x»,...^ can be written as the multidimensional probability density function for the vector x : p (a).

Independence Stochastic x

probability density function is given by:

Stochastic variables x.,x_,...,x are independent if the joint

P X l, x , x ( W - V ' 0 ^ = P x / V V ("2)""'Px ( V ( 2'3 7 )

1 2 n 1 2 n

Conditional probability density function

In many situations, actual information about the stochastic variables is available. Let x and y be stochastic vectors and the realization g of y be known. The conditional probability density function of x, given that the realization of y is 8, is defined by:

Px|v («|B> - P*,y (° 'B ) (2.38)

|y Py(B)

Throughout this thesis the probability density function p (fj) is assumed to be non-zero.

From Equation (2.38) Bayes theorem for probability density functions can be derived:

For two independent stochastic vectors it can be seen from Equations (2.37) and (2.38) that

Px|y ( a'8 ) = Px ( a ) ( 2 , 4 0 )

Stochastic process and stationarity

In geohydrology, many variables and vectors are functions of time. A stochastic variable or vector that is a function of time is called a stochastic process. The probability distribution and, consequently, the probability density function of a stochastic process are also functions of time.

A stochastic process is stationary if the probability distribution is independent of time.

Expectations

Mathematical manipulations with probability density function are often difficult. In most cases it is easier and sufficient to manipulate only with expectations of a stochastic variable, defined by:

-H»

E[ g(x) ] = ! g(x) p (a) da (2.41)

_ o o X

where: E[ . ] is the symbol for expectation

g(x) is a function of the stochastic variable.

The most important expectations are: - the expected value E[ x ]

- the variance E [ { X - E [ X ] Y ]

The concept of expectations can be also applied to multidimensional probability density functions and conditional probability density functions of stationary stochastic processes. For multidimensional stationary stochastic processes the expected value is a vector and the variance is extended to the covariance matrix:

E[ {x - E[ x ] }{ x - E[ x ] }T ]

The c o n d i t i o n a l expected v a l u e of x given y i s defined by:

E[ x|y ] = *2

P

| y

l

'

To solve practical problems it is often sufficient to do mathematical manipulations directly with expectations. The most important calculation rules for expectations are:

E[ c ] = c (2.43) E[ c x ] = c E[ x ] (2.44)

E[ x + y ] = E[ x ] + E[ y ] (2.45)

where x and y are stochastic variables or processes and c is a constant. If two stochastic variables x and y are independent, the rule

E[ xy ] = E[ x ] .E[ y ] (2.46) also is valid.

For example, the expected value of the state vector at time k in Equation (2.30) is

E[ 4>k ] = E[ A ^ + B uk + wk ] (2.47)

Because the elements of the matrices A and B are constants and the input vector u, is known, applying the calculation rules for expectations, it follows that:

E[ \ ] - A E[ *k_1 ] + Buk + E[ wk ] (2.48)

If the state vector at time (k-1) is known, the conditional expected value of the state vector at time k, given <(>. ., can be expressed by:

E

[

\ \ \ „ I ] "

*

_! +

+

[ "kK-l ^

'

vector at time (k-1) Equation (2.49) can be written as:

Multidimensional Gaussian distribution

The multidimensional Gaussian probability density function for the N-dimensional vector w, is defined as:

k

T

pw ( o )

N72 r

~

°-

[

] >• \

•

-

[

k 3

wk (2Tr)N/Z(det Q )! * k k

(2.51)

where: det Q, is the determinant of the covariance matrix Q, of the system vector w, .

Statistical analyses have been made of series of piezometric levels measured in 12 selected observations wells in The Netherlands (Van Geer 1981). The results of these analyses showed that the probability density functions of those measurement series are approximately Gaussian. So, in the present study, it is assumed that the state vector (piezometric levels) has a multimensional Gaussian distribution.

Three important properties of multidimensional Gaussian probability distributions are stressed here.

1. As can be seen from Equation (2.51) a multidimensional Gaussian probability distribution is completely determined by the expected value and the covariance matrix of the stochastic process.

2. Any linear transformation of a Gaussian vector results in a Gaussian vector.

3. The joint probability density function of two independent Gaussian processes is Gaussian.

From properties 2 and 3 it can be derived that if the initial condition <j>. and the system noise vector w. (i = l,2,...,k) have Gaussian distributions and are all independent, the state vector <j>. is also Gaussian distributed. By repeated use of Equation (2.30) for (k-1, k-2 1) the state vector at time k can be expressed as:

*k "

" *o

S J

k-i

V i >

-

(2.52) consists of k+1 terms, each of which containing one Gaussian vector (<j>n, w., w„,...,w,). Because of Property 2 these terms have

Gaussian distribution and from Property 3 it follows that the state vector (j>, also has a Gaussian distribution.

2.4 Measurement equation

The piezometric level is measured in a number of observation wells. Assume that the piezometric level is measured in observation well 1 at grid point (i,j). To relate the measurement to the state vector the latter is premultiplied by a row vector, which elements are equal to zero except for the element corresponding with grid point (i,j). This element is equal to one. Each measurement is corrupted with a measurement error. The relationship between measurement, the state vector and the measurement error at time k is given by:

jrl f k-[0...010...0]

n,i

♦I.J

+v

where: y. . is the measurement at Point (i,j) at time k [L] v} , is the measurement error at Point (i,j) at time k [L].

If there are m observation wells a measurement vector, (2.54), and a measurement error vector, (2.55), can be defined.

'l.k (2.54) 7

m,k J

The relationship given by Equation (2.53) can be generalized for m observation wells.

y

-V

v

(2.56)

where: C, is the measurement matrix. This matrix is determined by the locations of the observation wells.

Equation (2.56) is called the measurement equation. Although the locations of the observation wells will not change with time, the measurement matrix C, is a time-dependent matrix for two reasons:

1. If, at any moment of time, data are missing, the dimensions of Equation (2.56) can be reduced at that particular time by dropping the rows of y, , C, and v, of the equation that correspond with the missing data.

2. At any time the number of observation wells taken into account can be changed. This makes it possible to use measurement series with different frequencies in the same run (see Chapter 4 ) .

When the observation wells coincide with the grid points of the calculation grid, the elements of the measurement matrix C, can be derived easily. If the observation wells are located arbitrarily in the area, determination of the elements of the matrix C, is more difficult, but even in this case the linear relationship (2.56) holds

(see Appendix I ) .

The vector v, can be seen as a realization of a multidimensional k

stochastic process. If both the measurement equipment and the observer cause no systematic errors, the expected value of the measurement errors equals zero.

E[vJ = 0 (2.57)

The covarlance matrix of the measurement errors depends on the accuracy of the observations. This covarlance matrix can, therefore, be used as a criterion for the quality of the measurements. The flow diagram for a linear system with observations is given in Figure 2.9.

♦o

.1

'W-f

" '

*

Figure 2.9: Flow diagram for a linear system with observations.

2.5 RECAPITULATION

In this chapter it has been shown that the piezometric levels in a finite number of grid points in an area can be seen as the state vector of a linear system. The behaviour of the groundwater can be described as a process in such a system by:

the state equation: d>, = A <j> + B «k + w. (2.30)

the measurement equation: y = C, <f>. + v. (2.56)

The state equation and the measurement equation are the basic equations for the derivation of the filter algorithm which is discussed in Chapter 4.

Because the system noise and the measurement errors are realizations of stochastic processes, only an estimate of the state vector can be made. For this reason, the concept of expected values has been

introduced. In particular, the conditional expected value and the covariance matrix are of importance for the optimization of the groundwater monitoring network.

Chapter 3

THE OPTIMIZATION OF MONITORING NETWORKS

3.1 Introduction

Measurements of hydrologlcal variables are used to obtain information about how hydrologlcal processes vary in time and space. Measurements serve to estimate the values of hydrologlcal variables not only at the monitoring locations and times, but also at other locations and at other times. For example, data from gauges along a river are used to estimate the waterlevel everywhere along that river. In the last few decades, there has been an increasing interest in the evaluation and optimization of existing monitoring networks and in the optimal design of new networks. This is because of data management problems, which increase with increasing numbers of data, and to developments in systems analysis and computer facilities, which have provided operational optimization algorithms. Reviews of recent developments in hydrologlcal monitoring network design including extensive references are given in Moss (1982) and in Van der Made ed. (1986).

Most groundwater monitoring networks serve to gain insight into the behaviour of the groundwater system for many purposes, for example: domestic water supply, Irrigation and design of subsoil constructions (see Chapter 1 ) . Because piezometrlc levels reflect the state of the groundwater system (see Chapter 2 ) , the behaviour of the groundwater system can be deduced from piezometrlc levels. Therefore the design of a groundwater monitoring network can be restricted to the design of a network for monitoring piezometrlc levels. As stated in Chapter 1, the quality of a groundwater monitoring network is determined by the accuracy with which the piezometrlc levels at any point in space and time can be estimated from the measurements. With better estimates of piezometrlc levels a more efficient water management is possible. However, an increased quality of the monitoring network implies higher cost. In the ideal case the optimal groundwater monitoring network (density, frequency and equipment) follows from an optimization procedure which includes the cost of data collection and data management and the benefits (or losses) for all objectives. Section 3.2 discusses the design of such an ideal monitoring network.

Section 3.3 shows that the optimal groundwater monitoring network cannot be computed. Although optimization, as sketched in Section 3.2 is impossible in practice, the quality of a groundwater monitoring network can be analysed by adopting of a number of simplifications. The analysis gives guidelines for the monitoring network design. In Section 3.3 the simplifications used in the present study are discussed. The basic elements of the monitoring network analysis are the estimation of piezometric levels and the calculation of the accuracy of those estimates. Estimation of piezometric levels can be based on a number of estimation principles. The principles necessary for understanding the remainder of this thesis are defined in Section 3.4. A review of the most common methods in hydrological monitoring network design is given in Section 3.5. Finally in Section 3.6 the reasons for choosing Kalman filtering are presented.

3.2 Optimal design of groundwater monitoring networks

Insight into the behaviour of the groundwater is gained by estimating the state of the groundwater system (see Chapter 2 ) . In principle there are two different sources of information that can be used for the estimating piezometric levels (=state).

- measurements - models.

In groundwater systems, the measurements are observations of piezometric levels at a number of locations at a number of points in time. The model is the description of the behaviour of the groundwater. Because only mathematical models are under consideration here, the model is expressed in a vectorial equation (Chapter 2 ) . Measurements and model are used together in an estimation algorithm, with which estimates of the state (state vector) of the groundwater system can be calculated. In the estimation algorithm, the model provides the structure of the behaviour of the groundwater, while the measurements give information about the actual piezometric levels. The estimates differ from the real state because of estimation errors (see Figure 3.1), which are dependent on

the estimation algorithm the system noise

measurements measurement errors estimates estimation algorithm ■ estimation errors fnodet system noise

Figure 3.1 Measurements, model and estimation algorithm.

Because the real state is considered as a stochastic process (see Chapter 2) the estimation errors are also stochastic variables. The probability distribution of the estimation errors is dependent on the estimation algorithm and on the probability distributions of the system noise and the monitoring errors. Better estimates result in smaller variances of the estimation errors, thus smaller probabilities of large errors.

A very important assumption in groundwater monitoring network

optimization is, that: there exists a relationship between the

probability distribution of the estimation errors and the losses of the objective for which the estimates are used. A reduction of

estimation errors results in increasing benefits. Better estimates require more measurements, better locations or a better model. This means that, generally, a reduction of the estimation errors must be paid for by an increase in the estimation cost. So the optimal groundwater monitoring network would follow from the relationship between the estimation errors and the losses and the relationship between estimation errors and cost (see Figure 3.2).

cost tosses

. sum of losses

•losses

estimation error

If it is assumed that these relationships are known for all objectives. Then the optimal monitoring network would follow from an optimization process which optimizes simultaneously the network, the model and the estimation algorithm (see Figure 3.3).

losses monitoring network (density, frequency, equipment) probability distributions - system noise - measurement errors i optimization model cost - measurements - data management - calculations

I

Figure 3.3. Optimization process.

3.3 Groundwater monitoring network analysis

The optimization of a groundwater monitoring network can be performed if

a. the exact objectives,

b. the losses as a function of the estimation errors for all objectives,

c. the probability distribution functions of the system noise and the monitoring errors,

d. the cost of measurements, data management and calculations, are known.

In practice however, it is impossible to optimize a groundwater monitoring network in the way sketched in Figure 3.3. This is because:

1. The groundwater monitoring network is multi-purpose. In most situations monitoring series of several years or decades are

required. Therefore groundwater measurements made today will be used to support future studies. Hence not all loss functions can be known at the time the network is designed.

2. In a number of applications it is difficult or even impossible to quantify the losses. This is particularly so when environmental aspects are involved.

3. The probability distribution function of the system noise is always unknown.

4. Estimates of piezometric levels are calculated with an estimation algorithm, in which the measurements are combined with a model. Within certain margins a trade-off exists between the model chosen and the information obtained by measurements. This information is determined by the number of measurements and the locations of the observation wells. An advanced model with few measurements may result in the same estimation accuracy as a simple model with many measurements. Also the reverse can occur. Models with many parameters may require more data than simple

lumped parameter models for the same estimation accuracy. Consequently there may be a number of equivalent combinations of monitoring network and model.

5. Within certain margins, the monitoring network density and the measuring frequency can be traded off.

6. All costs and losses are time-dependent. The measurements are taken in the course of time. A model is developed at a particular point in time and the estimates are calculated, while the losses are related to the future. This makes it difficult to compare cost and losses at the time the monitoring network is designed.

7. In an area where nothing is known about the groundwater system, an initial groundwater monitoring network must be designed to gain insight into that system. After a period of time, measurements are available and the initial monitoring network can be evaluated, which may result in changes in locations of the observation wells or measurement frequencies.

Because of these seven reasons it is not possible to design the optimal groundwater monitoring network (Davis et. al. (1979), Dawdy

(1979), Moss (1979a), Tasker and Moss (1979). Instead, monitoring networks are designed by simplifying the optimization process. The ways in which the process can be simplified is dicussed below.

a. The unknown probability distributions of the system noise, the measurement errors and the estimation errors are characterized by expected values and covariance matrices of those variables. In most (or maybe all) studies this simplification is used. The advantage of considering these two quantities instead of the complete distribution is obvious. They can be estimated relatively easily from measurements (see Chapter 4 ) . Moreover, if the variables have Gaussian distributions, the expected values and covariance matrices contain all the relevant information. This is because a Gaussian distribution function is described completely by the expected values and the covariance matrix. The expected values and the covariance matrix can be interpreted physically. For an unbiased estimate of the state, the expected values of the estimation errors should be zero. Non-zero expected values of the estimation errors indicate an incorrect or poorly calibrated model. The covariance matrix (in particular the diagonal) is a measure for the scattering about the mean. If the probability distribution is (approximately) symmetrical, the covariance matrix is an appropriate criterion in monitoring network design. In the case of a Gaussian distribution function any confidence interval attached to the estimates can be calculated if the covariance matrix is known. If the probability distribution is asymmetrical, however, the estimates based on the expected value and the covariance matrix only should be interpreted very carefully.

Fortunately in groundwater flow, the probability distribution function of the state vector is often approximately Gaussian (Geer V. 1981). If the estimation algorithm used is linear, the estimation errors are also approximately Gaussian distributed (see Chapter 2) and hence the expected value and the covariance matrix contain most of the relevant information.