Identification of illegal groundwater pumping in semi-confined aquifers

Abstract

Illegal pumping from underground reservoirs is approached as an inverse problem. The procedure seeks among a set of suspected potential areas the ones from which water is effectively pumped out (or in) and at what rate does this take place? Influence coefficient and finite element are combined in order to simulate the aquifer response to hypothetical pumping scenarios. Then the optimal scenario which makes of simulated head distribution the best fit to the observed headset is retained as solution to the foregoing question. This optimal scenario is derived by minimizing the classical least squares error using the Levenberg-Marquardt Algorithm (LMA). Under the conditions of the treated aquifer example it is observed that LMA is non-convergent unless a regularizing measure is taken during the calculation of the correction applied to the parameter-vector between two successive iterations. The procedure is tested with uniform and non-uniform pumping. In both cases, it yields good identified pumping rates within the observation period, but its performance deteriorates markedly for the rates corresponding to the days when no monitoring well is yet active.

Citation Saffi, M. & Cheddadi, A. (2010) Identification of illegal groundwater pumping in semi-confined aquifers. Hydrol. Sci. J. 55(8), 1348–1356.

Résumé

L'identification des pompages illégaux à partir des réservoirs souterrains est approchée comme un problème inverse. La procédure consiste à chercher parmi une série de positions potentielles suspectes, celles d'où l'eau est effectivement pompée (ou injectée) et avec quel débit cela a-t-il lieu? Les coefficients d'influence et la méthode des éléments finis sont combinés pour simuler la réponse de l'aquifère à des scénarios de pompage hypothétique. Ensuite, le scénario optimal qui rend la distribution du potentiel hydraulique simulé la meilleure approximation de celle observée est retenu comme réponse à la question soulevée. Ce scénario optimal est obtenu indirectement en minimisant l'erreur au sens des moindres carrés en utilisant l'Algorithme de Levenberg-Marquardt (ALM). Sous les conditions de l'exemple considéré il est remarqué que l'ALM ne converge pas à moins qu'une mesure régularisante est prise lors du calcul de la correction appliquée au vecteur des paramètres entre deux itérations successives. La procédure est testée sur des scénarios de pompage uniforme et non uniforme. Dans les deux cas elle identifie correctement le débit pompé pour les jours d'observation, mais ses performances se détériorent nettement vis-à-vis des débits correspondants aux jours où les observations n'ont pas encore été commencées.

Key words:

• illegal pumping
• identification
• semi-confined aquifer
• inverse problem
• finite elements
• singular value decomposition (SVD)

Mots clefs:

• pompage illégal
• identification
• aquifère semi-confiné
• problème inverse
• éléments finis
• décomposition en valeurs singulières (DVS)

INTRODUCTION

Identification of illegal pumping, always carried out under camouflage, remains among essential measures for fair groundwater sharing. In practice, the authority managing an aquifer may have doubts about possible illegal groundwater abstractions. Typically, suspicion may grow around private property such as farms and factories where water pumping and wastewater injection is likely to occur. In such cases, a hydraulic head monitoring campaign can be run over a set of wells to help identify the locations from which the aquifer is probably being tapped.

From the governing equation viewpoint the recovery of the source term lies in the inverse problem category, and meets the second type of inverse problems in Sun's classification (Sun, 1999). The first and the third types are, respectively, identification of physical parameters such as transmissivity, and identification of initial and/or boundary conditions. Sun referred to the fourth type as any combination of the previous three.

It is noteworthy that identification of groundwater pollution sources has received much attention, and many aspects of the question have been debated (Gorelick et al., 1983; Aral & Guan, 1996; Neupauer et al., 2007). However, its counterpart, the problem of identifying illegal wells, still lags behind. This work attempts to contribute to addressing this issue by combining calibrated numerical models and observations. It extends the problem formulated and solved for a one-dimensional semi-confined aquifer by Saffi & Cheddadi (2007). So, given a reliable calibrated direct code and a record of observed heads, the point is to arrive at the well configuration (locations and magnitudes) which gives the best possible reproduction of the monitored hydraulic head.

The first section of this paper formulates the problem within an optimization framework. Then a test problem dealing with a two-dimensional (2-D) semi-confined aquifer is introduced, followed by a discussion. The latter is carried out first with uniform pumping, in order to illustrate through the considered example why the Levenberg-Marquardt Algorithm (LMA) diverges, and also to see which regularizing measure is required to make it converge. The discussion closes with the results of a non-uniform pumping and sensitivity analysis.

PROBLEM STATEMENT

Two kinds of wells are distinguished. Legal wells (LW), often used to satisfy municipal water demand, are known in number, location and magnitude. In the following they are represented by a vector , where Q _LWj (t) stands for the pumping rate of the jth legal well at time t, and is their total number.

Illegal wells are the areas from which water is abstracted or injected into the aquifer without prior permission. So, they are neither known in number nor in magnitude or location. However, here we refer to potential illegal wells (PIW), being the locations which are thought to be subject to illegal pumping. Similarly, PIWs are represented by a vector , where P _PIWj (t) denotes the magnitude of the jth PIW and is the number of such wells. The parameter-vector P (t) contains the problem unknowns to be determined.

Now, given a well field configuration ( P (t), Q (t)), the hydraulic head distribution h (t) over a flow region Ω is the solution to the following dynamical equation with the indicated initial condition (Pinder & Gray, 1977; Bear, 2007):

(1)

where h (t) = (h ₁(t), h ₂(t),…, h_n (t)) ^T ; h_j (t) is the hydraulic head in node j at time t; n is the number of discretization nodes; h ⁽⁰⁾ is the vector of initial heads; A and B are (n × n)-matrices written in terms of the aquifer storativity and transmissivity, respectively; and F is the forcing vector containing point and distributed sources/sinks and the boundary conditions. The values of A , B and F are computed for the 2-D semi-confined flow equation recalled in Appendix A, and discretized using conventional linear triangular finite elements. For convenience, equation (1)(B1) is presented as an input–output relationship highlighting the parameter-vector P (Saffi, 2008):

(2)

In equation (2)(B2), h ^(s) = (h ₁ ^(s), h ₂ ^(s), …, h_n ^(s)) ^T , where h_j ^(s) denotes the hydraulic head in node j and at time s; is the vector of nodal hydraulic head when all illegal wells are inactive; Δt is the discretization time step; and G ^(s) stands for the matrix of influence coefficients at time s. The vector contains the known pumping rates with Q _LWj ^(k) referring to the jth LW at time k and s _max being the number of time levels over which the system is considered. Likewise, the vector represents the unknown pumping rates where P _PIWj ^(k) refers to the jth PIW during time k. Practically, G ^(s) is computed only once by solving system (1), in which boundary and initial conditions are zeroed; and all the sources and sinks are replaced with a unit Dirac delta function. Then, G ^(s) is stored in a separate file ready-for-use alongside equation (2)(B2) for any potential scenario P . In this way, a great deal of computation time is saved during inverse problem resolution. A full account of this is found in Illangasekare & Morel-Seytoux (1982). Explanatory details on the structure of the matrix G ^(s) are also presented in Appendix A.

The aquifer system is observed through monitoring wells MW1, MW2,…, MW, over a period of time [s _min, s _max]; where s _min corresponds to that time step when the first head observation is made. Let [h _MWj ^obs]^(s) denote the observed head at monitoring well MWj and time step s for and s _min ≤ s ≤ s _max. For convenience, the [h _MWj ^obs]^(s) terms are gathered in the vector defined by:

(3)

In parallel, the simulated hydraulic head at MWj and time step s, i.e. the component h _MWj ^(s) of the vector h ^(s) of equation (2)(B2), will be denoted [h _MWj ^sim]^(s). Again, the [h _MWj ^sim]^(s) terms are put in the following vector:

(4)

Next, a feasible model input P which achieves the best fit between the foregoing vectors h ^sim( P ) and h ^obs, is sought by minimizing the criterion:

(5)

TEST PROBLEM DEFINITION

The following example is inspired by Willis & Yeh (1987). Figure 1(a) displays a hypothetical flow region Ω delimited by two river reaches that impose a constant head of 100 m. The two other (shaded) sections of Ω border are of no-flow type. The aquifer transmissivity and storage coefficient are 1000 m² d^‐1 and 0.01 respectively. Semi-confining layer hydraulic conductivity and thickness, K_a and m_a , are set equal to 0.01 m d^‐1 and 30 m, respectively. The external hydraulic head H_a  = 100 m. The flow governing equation is recalled in Appendix A. Water is pumped from the only legal well LW1 located in the node 28 at a daily rate of 10 000 m³ (i.e. , the number of legal wells is equal to 1). Steady-state conditions are assumed to prevail then.

Fig. 1 (a) Aquifer system configuration: black circle denotes the only legal well (LW1); stars indicate potential illegal wells (PIW); unfilled circles show monitoring well (MW) positions. (b) Initial hydraulic head distribution taken as the steady-state solution to equation (1)(B1) for: Q _LW1 = 10 000 m³ d^‐1; and P _PIW1 = P _PIW2 = 0.

At some time level s = 0, two areas around nodes 34 and 44 are suspected to be subject to illegal pumping (i.e. , the number of potential illegal wells is 2). That is, a monitoring campaign was started in monitoring wells MW1, MW2,…, MW6 located in nodes 26, 36, 43, 25, 50 and 6, respectively, over a 9-day-period of observation that starts at s _min = 7 and ends at s _max = 15 With these data there are more equations than unknowns and problem (5) is over-determined. Synthetic observed heads are obtained by perturbing to 5% the solution of equation (2)(B2) started from the initial head distribution of Fig. 1(b), with the flow domain segmented into 97 triangular elements using n = 62 nodes and a time step Δt = 1 day.

RESULTS AND DISCUSSION

In the following, the Levenberg-Marquardt Algorithm (LMA) is used to solve minimization problem (5) for the above-defined test problem; first with a uniform illegal pumping of 2500 m³ d^‐1 in magnitude for both PIW1 and PIW2 and then with a non-uniform pumping scenario.

After launching the algorithm from several starting points, it is observed that J is nearly constant and takes very small values while Δ P ^(k) remains large. The term Δ P ^(k) is the correction applied to the parameter-vector P ^(k–1) during the kth iteration; it is calculated via equation (B1) (see Appendix B). For instance, the last row of Table 1 reports that J is as small as 3 × 10^‐6 while Δ P ^(k) is of the order 13 500 m³ d^‐1: the algorithm diverges. The cause underlying this poor performance lies in the relationship between Δ P ^(k) and the subsequent variation J( P ^(k)) – J( P ^(k–1)) of criterion J denoted ΔJ ^(k). From the Taylor series applied to J( P ^(k) =  P ^(k–1) + Δ P ^(k)), it may be shown that:

Table 1  The first five iterations of non-regularized LMA for PIW1. Similar results (not reported here) are obtained for PIW2

P ^(0)	Δ P ^(1)	P ^(1)	Δ P ^(2)	P ^(2)	Δ P ^(3)	P ^(3)	Δ P ^(4)	P ^(4)	Δ P ^(5)
0	−2.3903	−2.3903	+0.0284	−2.3619	+0.1003	−2.2616	+0.5408	−1.7208	+2.9714
0	−2.4365	−2.4365	−0.0436	−2.4801	−0.1299	−2.6099	−1.0529	−3.6628	−6.0751
0	−2.4914	−2.4914	−0.0546	−2.5460	+0.0578	−2.4882	+0.0853	−2.4029	+0.8793
0	−2.5349	−2.5349	−0.0539	−2.5888	+0.1832	−2.4056	+0.7486	−1.6570	+3.8010
0	−2.5198	−2.5198	−0.0523	−2.5721	+0.0159	−2.5563	−0.1881	−2.7444	−2.1994
0	−2.4076	−2.4076	+0.0030	−2.4046	−0.1110	−2.5156	−0.1013	−2.6169	+0.3898
0	−2.5397	−2.5397	+0.0386	−2.5011	+0.0289	−2.4722	+0.0223	−2.4499	−0.0342
0	−2.4834	−2.4834	−0.0247	−2.5081	−0.0123	−2.5204	−0.0087	−2.5291	−0.0015
0	−2.4938	−2.4938	+0.0087	−2.4852	+0.0042	−2.4809	+0.0019	−2.4790	−0.0005
0	−2.5060	−2.5060	−0.0060	−2.5119	−0.0021	−2.5140	−0.0014	−2.5155	−0.0017
0	−2.4808	−2.4808	+0.0008	−2.4800	+0.0009	−2.4790	+0.0003	−2.4788	−0.0003
0	−2.5005	−2.5005	−0.0022	−2.5027	−0.0003	−2.5030	−0.0003	−2.5033	−0.0008
0	−2.4973	−2.4973	−0.0007	−2.4980	+0.0003	−2.4977	+0.0001	−2.4976	−0.0004
0	−2.5034	−2.5034	−0.0017	−2.5051	+0.0001	−2.5050	−0.0000	−2.5051	−0.0005
0	−2.4907	−2.4907	−0.0002	−2.4909	+0.0001	−2.4908	+0.0000	−2.4908	−0.0003
14.48^(a)	13.78^(b)	5.08 × 10^‐6 (a)	0.485^(b)	3.09 × 10^‐6 (a)	1.02^(b)	2.99 × 10^‐6 (a)	3.15^(b)	2.92 × 10^‐6 (a)	13.66 ^(b)
^(a)Value of criterion J in P ^(k).
^(b)Δ P ^(k) Euclidean norm magnitude (in 10^3 m^3 d^‐1).
Entries of vector P ^(k) are in 10^3 m^3 d^‐1.

(6)

where the dot denotes the scalar product; and a ^(k) is a vector the exact expression of which is not relevant in the present context. From equation (6)(6) it is clear that the part Δ P ₂ ^(k) of Δ P ^(k), which is nearly perpendicular to a ^(k) brings no substantial contribution into ΔJ ^(k). In other words, the variation in J stems essentially from the difference (Δ P ^(k) – Δ P ₂ ^(k)) termed Δ P ₁ ^(k). So to gain more insight into why J stays almost constant, magnitudes of Δ P ₁ ^(k) and Δ P ₂ ^(k) were closely watched by running the code from several starting points. It is observed that vector Δ P ₂ ^(k) (i.e. the one with which criterion J does not vary noticeably) is behind the large length of the correction-vector Δ P ^(k). This is clearly conveyed by the Euclidean norm displayed on the second row of Table 2. In contrast, the useful component Δ P ₁ ^(k) which instigates the variation of J shrinks dramatically with iterations (first row of Table 2). Consequently, J stands nearly constant, while Δ P ^(k) grows out of proportion. These conflicting roles played by Δ P ₁ ^(k) and Δ P ₂ ^(k) make the algorithm definitely non-convergent. So a practical step to take while playing the LMA is to evaluate Δ P ₂ ^(k) and clear it from Δ P ^(k) during calculation. This regularizing measure is tested in the following numerical experiments.

Table 2  Comparison of Δ₁ ^{(k )} and Δ₂ ^{(k )} Euclidean norms (in 10³ m³ d^{‐ 1}) for the five iterations displayed on Table 1

	Iteration 1	Iteration 2	Iteration 3	Iteration 4	Iteration 5
||Δ P 1 ^(k)||	13.7131	0.4700	0.0514	0.0041	4.95 × 10^‐4
||Δ P 2 ^(k)||	1.3240	0.1178	1.0155	3.1504	13.6622
||Δ P ^(k)||	13.7769	0.4845	1.0168	3.1504	13.6622

Now, to turn the previous conceptual reasoning into computer code, an explicit means to quantify Δ P ₂ ^(k) is required. This means is provided by the singular values spectrum of the curvature matrix C of equation (B1). Paragraph 2.3 of Nazareth (1980) adduces mathematical grounds for this. The complete orthogonal factorization of the Jacobian matrix used then should be replaced with the following singular-value decomposition of the matrix C :

(7)

where U is a column-orthogonal matrix of the same dimension as C ; W is a square diagonal matrix with nonnegative elements w_i referred to as C ’s singular values and V is a square column-orthogonal matrix. Practically, matrices U , V and W are returned by the standard MATLAB command svd(), or Singular Values[] for MATHEMATICA users. Package svdcmp (for PASCAL, FORTRAN and C languages) doing a similar job is provided by Press et al. (1989).

In all, up to some cut-off order M ₀ the w_i values stay with a same order of magnitude, but beyond the cut-off they drop dramatically. A typical profile of matrix C singular values is depicted in Fig. 2 with a corresponding M ₀ of 23. In this case C comes from the aquifer example of Fig. 1 under uniform pumping. Regression lines A and B in Fig. 2 with respective slopes of −0.13 and −1.25 illustrate how sharp is the decline in the magnitude of w_i ; M ₀ is the threshold which demarcates Δ P ₁ ^(k) from Δ P ₂ ^(k) in the following way:

(8)

Fig. 2 Typical profile of singular values w_i of the curvature matrix C resulting from the aquifer example of Fig. 1 under uniform pumping of 2500 m³d^‐1and a data set perturbed to 5%. The regression lines A and B with slopes −0.13 and −1.25 illustrate how important is the decline in w_i magnitude. The black arrow points at the first substantial drop in two successive w_i , i.e., w ₂₃ and w ₂₄.

The foregoing Δ P ₁ ^(k) expression meets equation (14.3.17) in Press et al. (1989). Vector b originates in the classical equations underlying the LMA. It is made explicit in Appendix B. The terms U _*j and V _*j are the jth column vectors of matrices U and V , respectively. Their numerical values are known from the above mentioned commands and procedures.

At this stage, a definite choice of M ₀ remains in order to make equation (8)(8) fully operational. Press et al. (1989) pointed out that no standard rule is available for this task, and the user is called to exercise some of his discretion. Actually the value of M ₀ is ad hoc to the treated technical problem. Practically, through the aquifer example under work, a rough range of three or at most four potential values of M ₀ is unmistakable; afterwards the exact value can be efficiently determined by trial and error. Yet, even a naked-eye study of Fig. 2 shows that the first large gap between two successive values occurs at w ₂₃.

Actually, () corresponds to the number of equations in system (B1) which are nearly redundant, so they are at best useless during Δ P ^(k) calculation. They are so because of the information lost before the head monitoring campaign was started, i.e. (s < s _min). This is again confirmed by heuristic code runs conducted with various s _min values. It seems that the smaller is s _min the larger is M ₀. A good preliminary estimate of M ₀ for this problem would be the dimension of the parameter-vector (), less the number of days during which no head observation was made (i.e. s _min).

Numerical experiments under regularized Δ P ^(k)

Now the LMA is rerun taking care at each iteration to drop Δ P ₂ ^(k) while solving system (B1). The new estimate P ^(k–1) of the parameter-vector and the corresponding regularized correction Δ P ^(k) are listed in Table 3. Notice that both J( P ^(k–1)) and ||Δ P ^(k)||₂ decrease correspondingly, and overall P ^(k–1) converges to the correct rate of 2500 m³ d^-1.

Table 3  The first five iterations obtained with regularized LMA (M ₀ = 23) for PIW1

P ^(0)	Δ P ^(1)	P ^(1)	Δ P ^(2)	P ^(2)	Δ P ^(3)	P ^(3)	10^4× Δ P ^(4)	P ^(4)	10^4× Δ P ^(5)
0	−2.4725	−2.4725	+0.0196	−2.4528	−0.0009	−2.4537	+0.4397	−2.4537	+0.0614
0	−2.5893	−2.5893	+0.0193	−2.5700	−0.0005	−2.5705	+0.2944	−2.5705	+0.0400
0	−2.6597	−2.6597	+0.0134	−2.6463	−0.0003	−2.6467	+0.1543	−2.6466	+0.0215
0	−2.6343	−2.6343	−0.0016	−2.6359	−0.0006	−2.6365	+0.0425	−2.6365	+0.0104
0	−2.4300	−2.4300	−0.0266	−2.4566	−0.0012	−2.4578	−0.0348	−2.4578	+0.0074
0	−2.4046	−2.4046	−0.0422	−2.4468	−0.0012	−2.4480	−0.1517	−2.4480	−0.0065
0	−2.4980	−2.4980	+0.0229	−2.4751	+0.0023	−2.4728	+0.0226	−2.4728	−0.0200
0	−2.5013	−2.5013	−0.0095	−2.5109	−0.0011	−2.5120	+0.0189	−2.5120	+0.0138
0	−2.4838	−2.4838	−0.0014	−2.4852	+0.0005	−2.4847	−0.0282	−2.4847	−0.0086
0	−2.5091	−2.5091	−0.0017	−2.5108	−0.0002	−2.5111	+0.0101	−2.5110	+0.0035
0	−2.4788	−2.4788	−0.0017	−2.4805	+0.0001	−2.4804	−0.0100	−2.4804	−0.0022
0	−2.5011	−2.5011	−0.0014	−2.5025	−0.0000	−2.5025	+0.0011	−2.5025	+0.0006
0	−2.4969	−2.4969	−0.0013	−2.4982	+0.0000	−2.4982	−0.0035	−2.4982	−0.0006
0	−2.5036	−2.5036	−0.0016	−2.5051	−0.0000	−2.5051	−0.0009	−2.5051	+0.0000
0	−2.4906	−2.4906	−0.0004	−2.4910	+0.0000	−2.4910	−0.0014	−2.4910	−0.0002
14.48^(a)	13.71^(b)	4.40 ×10^‐6 (a)	0.15^(b)	3.41 ×10^‐6 (a)	0.010^(b)	3.41 ×10^‐6 (a)	1.5 ×10^‐4 (b)	3.41 ×10^‐6 (a)	2.72 ×10^‐5 (b)
^(a)Value of criterion J in P ^(k).
^(b)Δ P ^(k) Euclidean norm magnitude (in 10^3 m^3 d^‐1).
Entries of vector P ^(k) are in 10^3 m^3 d^‐1.

Figure 3 depicts the identified pumping rates for both PIW1 and PIW2 over the considered 15-day period. It appears that the pumping rates corresponding to days s, between s _min = 7 and s _max = 15, are well identified for both wells, while those relative to days s out of the observation period (s < s _min) are poorly identified. Moreover, because of the storage effect (aquifer memory), observed heads during early days of the campaign (e.g. s = 7 or 8) still carry some information on the pumping rates of the latest unmonitored days (e.g. s = 5 or 6). So, depending on how long/short the system memory is, the general trend is this: prior to s _min the farther one goes back in time the less accurate the identified pumping becomes, as it is exemplified on PIW2 in particular during the first three or four days (Fig. 3).

Fig. 3 Identified pumping rates for the aquifer of Fig. 1 under a uniform pumping of 2500 m³ d^‐1. Data are perturbed to 5%. Large oscillations are shown on the first 6 days when the hydraulic head is not monitored yet.

Figure 4 reports the results of a second example: data are perturbed to 5% and the methodology is similar to the previous case except that the pumping magnitudes at PIW1 and PIW2 are made to vary with time as shown by the dashed line. Again LMA performs well within the observation period [s _min, s _max] but it poorly identifies the pumping rates lying outside of it. Moreover, a similar test to that carried out by Saffi & Cheddadi (2007) in which data are perturbed to 2.5%, 5% and 7.5%, shows that the regularized pumping rates are not obscured by the added noise though the quality of the results slightly deteriorates.

Fig. 4 Identified pumping rates for the aquifer of Fig. 1 under a non-uniform pumping. Data are perturbed to 5%.

CONCLUSION

This work formulates and solves the problem of illegal pumping in 2-D semi-confined aquifers. It aims to determine pumping rates and the locations from which groundwater is being illegally abstracted among a known and finite set of suspected areas.

The question comes down to identifying the vector P of hypothetical pumping rates P _PIWi ^(s) at potential illegal well number i and time step s, 1 ≤ s ≤ s _max, from a set of observed heads [h _MWj ^obs]^(k) carried out at monitoring well j and time step k (s _min ≤ k ≤ s _max).

Parameter-vector P is determined by minimizing the criterion J of equation (5)(5) using the LMA. As a first approach the procedure is explained using a hypothetical aquifer example in which it is assumed that the main parameters of the system alongside initial and boundary conditions are accurately known. Under these conditions it appears that:

1. The problem is well-posed with respect to magnitudes P _PIWi ^(s) corresponding to day s for the monitoring campaign (s _min ≤ s ≤ s _max), while it is ill-posed for P _PIWi ^(s) with s < s _min.

2. The LMA converges if, for each iteration, only the most significant M ₀ singular values w_i of the curvature matrix are taken into account during the calculation of the correction applied to the parameter-vector.

3. A rough surrounding of the cut-off order M ₀ can be delineated from the w_i profile; then the exact value is practically determined by trial and error. It is observed that M ₀ is strongly linked to s _min, the number of days when no head observation was made. The shorter s _min, the larger is M ₀ and the more well-posed is the problem.

However, the problem should be investigated under less restricting hypotheses and the effect of error propagation in the estimation process needs assessment.

APPENDIX A

Semi-confined flow governing equation

The mathematical model considered in this paper rules 2-D semi-confined flows. It is written (Bear, 2007):

(A1)

where h is the hydraulic head, T is the aquifer transmissivity; K_a and m_a are respectively the conductivity and the thickness of the layer overlying the aquifer; H_a is the external hydraulic head; and is the storage coefficient.

Matrix of Influence coefficients, G ^(s)

In this work, the matrix G ^(s) of equation (2)(B2) comes from finite element discretization of the above equation (A1). It is a concatenation of sub-matrices ( is the number of PIWs). Each sub-matrix is associated with a PIW and contains s _max columns and n rows, as explained by equation (A2); s _max is the number of time levels over which the aquifer system is considered; and n is the number of grid nodes.

(A2)

Consequently, G ^(s) is a array. The general term G _p,PIWk ^s,i is read as the aquifer response at time level s, in node p to a unit stress applied at time level i in the PIWk under homogeneous boundary and initial conditions. In other words, G _p,PIWk ^s,i is the following derivative of the hydraulic head with respect to pumping:

(A3)

Matrix G ^(s) is relative to the present time step s, so its columns beyond (s + 1) are conventionally nil (i.e. G _p,PIWk ^s,i  = 0 for i ≥ s + 1). A detailed procedure to compute the influence coefficients is found in Maddock (1972).

APPENDIX B

Computation of correction Δ P

The correction Δ P applied to parameter-vector P while minimizing the criterion J of equation (5)(5) using the Levenberg-Marquardt algorithm is solution to the following set of equations (Press et al., 1989):

(B1)

The curvature matrix and the vector are defined for the criterion J by:

(B2)

(B3)

P_i is a provisional notation of the ith component of the parameter-vector P as introduced in equation (2)(B2); δ _ij is the Kronecker symbol; and λ is an adaptable parameter. Further details are given by Press et al., 1989, pp. 574–576).