4. Eigenmodels

The conference paper in Appendix I, about the eigenvalue approach to modelling aquifers, provides the theoretical basis for this section.

4.1 Assumptions

Figure 2 illustrates the concepts and mathematical symbols used in the following discussion.

Figure 2: Variables and parameters for eigenmodel assumptions

4.1.1 Dynamics of Recharge Processes

There is a fundamental difference in the responses of piezometric levels in an aquifer to changes in river recharge and land surface recharge.

Time-varying piezometric effects near a river, that are associated with changes in river recharge, are rapidly attenuated with distance from the river into the aquifer. The time lag for piezometric effect also increases with distance from the river.
Changes in land surface recharge cause relatively rapid piezometric effects everywhere in the aquifer, and the magnitude of these effects increases with distance from the fixed-head boundaries of the aquifer.

Therefore we assume that at any observation well, the piezometric effect of river recharge r(x,y) is a constant value for that location and that temporal variations in piezometric head u(x,y,t) are caused by the effects of land surface recharge and time-varying abstraction. The validity of this assumption improves with distance from river recharge sources.

4.1.2 Land Surface Recharge Pattern

Land surface recharge is assumed to have a fixed spatial distribution, but the overall magnitude varies with time. This is expressed mathematically by stating that the recharge RLS(x,y,t) at time t for location x, y is given by:

(9)

where P(x,y) is a spatial distribution pattern and R(t) is a time-varying magnitude. This means that the output from one water balance model may be quite satisfactory for a region with varying amounts of rainfall because the spatial variations are incorporated into the eigenmodel calibration, as will be shown later.

4.1.3 Aquifer Properties

There are no assumptions required about the variation of transmissivity T(x,y) and storativity S(x,y) throughout an aquifer. The geological structure may be quite heterogeneous. However, the assumption of 2-D groundwater flow means that piezometric levels in areas with a significant vertical component of groundwater flow may be less well estimated.

4.1.4 Dynamic Effect of the Vadose Zone

Land surface recharge must travel through the unsaturated region (vadose zone) between the land surface and the groundwater before it has any piezometric effect. Recharge into the top of the vadose zone displaces water from the bottom of the zone into the groundwater surface, by means of hydraulic wave propagation. This process introduces a time delay and some attenuation of the estimated land surface recharge signal. We include an additional linear storage in the eigenmodel to account for this effect. This storage element can also account for the dynamic effect of groundwater perching above an aquitard as it leaks into a semi-confined aquifer.

4.2 Model Structure

An eigenmodel of the dynamic behaviour of an aquifer can be represented as a set of linear storages connected in parallel, as shown in Figure 3. The dynamic effect of the vadose zone is represented by a single linear storage in series with the parallel set.

Figure 3: Eigenmodel structure

4.2.1Characteristics of an Eigenmodel

The eigenmodel structure shown in Figure 3 has some important characteristics that are relevant to practical applications:

The dynamic parameters k1, k2 ..., which are the eigenvalues of this linear system, are the same at all locations in an aquifer. This property enables transfer of these parameter values from locations with good data to other locations with poorer data.
Although the number of parallel elements is theoretically large, in practice only a few elements are required. Our experience is that for data at monthly time intervals only three need be considered.
The smallest eigenvalue k1 has the most influence on long-term dynamic behaviour, and the other, larger, eigenvalues account for the more rapid response of piezometric head to changes in recharge.
The larger eigenvalues are more likely to be significant near fixed-head boundaries, because they can account for the rapid dynamics of the shorter drainage paths in these regions.
The vadose zone dynamic parameter kv does vary with location, but its effect is usually less significant than the smallest eigenvalue, so it can be neglected (set to an arbitrary large value) during the initial stage of model calibration.
The gain coefficients gi(x,y) do depend on location.

4.3 Equations for Spreadsheets

4.3.1 Prediction Equation

Differential equations, similar to equation (1), can be written for each of the components of the eigenmodel structure shown in Figure 3 but these would not be suitable for direct solution in a spreadsheet. We need a difference equation, similar to equation (6), which can be used with discrete time-series data such as monthly totals of recharge and monthly values of groundwater levels. Equation (6) is a first-order difference equation, because it has only one "a" coefficient. The four-component structure of Figure 3 will yield a fourth-order difference equation by use of a mathematical technique called the method of transforms (in this case, z-transforms). The development of the equation is shown in Appendix II, and the resulting form is:

(10)

where ĥ_n is the prediction of the observed piezometric head h_n.

4.3.2 Parameter Relationships

The parameters ai and bi of equation (10) do not have obvious physical meanings, because they are derived from the more physically relevant parameters ki and gi(x,y) by the mathematical transformation processes. The parameters in equation (10) are also interdependent, which can be a problem for model calibration, and it is difficult to set meaningful initial trial values during calibration. Therefore, calibration is conducted on the z-transform model, which has parameters αi and βi, and these values are converted to the ai and bi of the difference equation by the following relationships (developed in Appendix II) in the spreadsheet:

(11)

(12)

The calibrated values of αi and βi are also converted to the equivalent eigenmodel parameters ki and gi(x,y) by the relationships:

(13)

(14)

Equations (11) to (14) are shown above only for the purpose of illustrating the analytical process, and are already embedded in the example spreadsheets. Figure 4 shows the results header from an example spreadsheet.

Figure 4: Results header from an eigenmodel spreadsheet example

The parameter values for the transform model, in the shaded line, have been obtained by calibration and the other parameter values have been calculated by means of equations (11) to (14).

4.4 Model Calibration

Eigenmodels expressed in spreadsheet form are calibrated by use of the Solver routine in the Tools menu of Microsoft Excel. There are no problems with having missing data in the time-series of observed piezometric head, other than some loss of information, because the objective function in the spreadsheet ignores blanks.

4.4.1 Objective Function

The Tools user-form prompts for a Target Cell, and this should be set to the cell Obj. fn. shown by dark shading in Figure 4, and the option Min selected. The objective function is the sum of squares of the model error. This error is called noise Nn and is defined as the difference between the observed and predicted values of piezometric levels:

(15)

Parameter Constraints

In the Subject to constraints box of the user-form, the Transform Model parameters (shaded line) should be set as follows:

Alpha-parameters: >=0, <=1
Beta-parameters: >=0
Base: no constraint

4.4.3 Initial Parameter Values

All the Transform Model parameters should be set to have an initial value of zero. This convenient setting is one of the advantages of conducting calibration on the z-transform version of the eigenmodel.

4.4.4 Calibration Procedure

The difference equation parameters are arranged in a particular order in the spreadsheet to facilitate a stepwise calibration procedure that minimises ambiguity and instability of parameter values. The Solver tool should be applied successively to the following selections of parameters shown in the shaded line of Figure 4:

Base
Base, Alpha1, Beta1
Base, Alpha1, Beta1, Alpha2, Beta2
Base, Alpha1, Beta1, Alpha2, Beta2, Alpha3, Beta3
The best of the above selections plus Alpha4

The above groupings are based on selecting, in a stepwise manner, parameters in order of their likely effect on variations in piezometric levels. The inclusion of Alpha4 is to apply the vadose zone element to the best groundwater model.

Contact for Enquiries

MAF Information Services
Pastoral House
25 The Terrace
PO Box 2526
Wellington, NEW ZEALAND

Fax: +64 4 894 0721
Contact this person