Use of the Kalman Filter for the Interpretation of Aquifer Tests Including Model and Measurement Errors

: The hydraulic parameters representative of actual aquifer conditions can be obtained through aquifer tests formerly known as pumping tests. Diverse methodologies based on analytical or numerical solutions have been proposed for the interpretation of aquifer tests; however, measurement and model errors are often neglected, which could lead to hydraulic parameter values that do not reﬂect the aquifer conditions. In this paper, a new alternative is presented for the interpretation of aquifer tests in conﬁned aquifers based on the Cooper–Jacob solution by means of the dynamic Kalman ﬁlter and a nonlinear optimization method. This proposal was tested in two previously published case studies; the measured drawdowns were ﬁltered by considering measurement and model errors to match the Cooper–Jacob solution. For the case studies, the results show that ﬁltering the measured drawdowns leads to variations of up to 49.97% in the values for T and 150% for S when compared to the values determined by methodologies that neglect measurement and model errors. A poor match between ﬁltered and measured data reﬂects large measurement errors and considerable deviations of the aquifer conditions with respect to the proposed model.


Introduction
It has been stated that numerical flow models represent the most advisable and powerful available tool for the adequate evaluation, planning and management of groundwater [1]. Most groundwater models are based on mass conservation and Newton's second and third laws (momentum and energy), and their analysis requires the solution of differential equations. To represent the heterogeneity of aquifer properties, the spatial distribution of various hydraulic parameters needs to be established, such as: the hydraulic conductivity (K, L 1 T −1 ), the hydraulic transmissivity (T, L 2 T −1 ), particularly for horizontal flow, the storage coefficient (S, L 3 L −3 ), the specific storage (Ss, L −1 ), and the hydraulic diffusivity (D, L 2 T −1 ), among others [2].
In order to calibrate a numerical flow model, the hydraulic parameters can be obtained from aquifer tests, formerly known as pumping tests, which evaluate the response of an aquifer (time-drawdown data usually measured at observation wells) by extracting water (at a constant flow rate) through a pumping well. Another option to determine the hydraulic parameters is the execution of permeability tests performed in the laboratory, or alternatively with the use of tracers, from reference tables, or empirical formulas; however, aquifer tests provide the most representative values for specific aquifer conditions [3][4][5].
A number of approaches have been proposed for the interpretation of aquifer tests; traditional methods based on analytical solutions for the interpretation of aquifer tests consist of looking for the best match between a type-curve and plotted drawdown measurements to obtain values that are substituted into simplified mathematical forms of the solution to determine the hydraulic parameters. The type-curve for a specific flow solution represents a function of the theoretical relationships among the hydraulic parameters, the elapsed time from the start of pumping and the distance between the pumping and the observation well. Some classical analytical solutions for transient flow have been developed for: (a) confined aquifers [6,7]; (b) unconfined aquifers [8,9]; (c) leaky confined aquifers [10][11][12]; and (d) fractured aquifers [13,14]. They represent the exact solution of a mathematical model that describes the groundwater flow in the vicinity of a pumping well following particular assumptions.
A variety of commercial software solutions based on analytical solutions have been developed, such as Aquifer win32 [15], AQUIFERTEST [16], AQTESOLV [17], ANSDIMAT [18], and SATEM 2002 [19]; typically, field data are matched to type-curves through an optimization method. Poor adjustment can be related to the fact that the actual conditions of the pumping well and the aquifer do not satisfy the assumptions considered to obtain the analytical solution; however, in other cases, this may be associated with measurement errors.
The approach is different when using numerical models for the interpretation of aquifer tests; a model is constructed and calibrated by evaluating different values for the hydraulic parameters to produce drawdown estimates that approximate the measurements [20,21]. Even though numerical models represent a more general alternative for the interpretation of aquifer tests due to their ability to consider additional field conditions that analytical solutions neglect, the methodologies based on the latter continue to be the most practical alternative for the determination of the hydraulic parameters if field conditions do not deviate significantly from the assumptions employed for the derivation of the respective solution.
In both approaches (analytical and numerical), measured drawdowns differ from the estimated values, using the determined hydraulic parameters, mainly due to measurement errors and deviations from the model assumptions; in this sense, these differences can be reduced by filtering noisy measurements through data assimilation which consists of combining a set of discrete measurements and a prediction model to obtain the best estimates for a state variable.
The Kalman filter has been used for data assimilation in hydrogeology in different applications, such as: the calibration of groundwater numerical flow models [22], the optimal design of monitoring networks [23][24][25], the estimation of hydrogeological variables [26,27], and the interpretation of aquifer tests [28,29]. Recent studies proposed the use of artificial neural networks to estimate hydrogeological variables [30].
Sen [28] applied the Kalman filter for the interpretation of aquifer tests in confined aquifers using the Theis equation; this updates the estimation of T and S (and their respective estimate error variances) each time newly recorded time-drawdown data are input within the procedure. Leng et al. [29] used the extended Kalman filter (EKF) to interpret aquifer tests for confined and unconfined aquifers; a cubic spline was used to estimate values at regular time intervals, which were used to estimate the hydraulic parameters. An important advantage of both procedures is that the length of time for aquifer tests can be reduced since the employed adaptative processes help to determine the hydraulic parameters using only part of the observed drawdown data; however, as with traditional methods for the interpretation of aquifer tests, the drawdown measurement errors are neglected, which in some cases could lead to high uncertainty in the results, providing values of the hydraulic parameters that would not be representative of the aquifer in the vicinity of the pumping well where the aquifer test was carried out.
In this paper, the Cooper-Jacob solution [7] was selected as the mathematical flow model to implement the dynamic Kalman filter [31]. This model for confined aquifers express drawdowns as a function of the elapsed time from the start of pumping; the constant values for this relationship are the flow rate in the pumping well, the distance between the pumping and the observation well, and the hydraulic parameters (T and S). Theoretically, the drawdown estimates obtained when filtering the observed data (considering the correct measurement and model errors) would produce values similar to the drawdowns obtained with the Cooper-Jacob solution. This procedure was used to calibrate the model (determine T and S) and to evaluate uncertainties during the execution of aquifer tests; the results were compared to those obtained following different approaches for two previously published case studies.

Materials and Methods
The Kalman filter is a sequential mathematical procedure for data assimilation that operates through a prediction and correction mechanism. The filter is sequential because it recalculates the solution each time a new measurement is available without using old data again. This procedure obtains a new estimate of the state from its previous estimate by adding a correction term that incorporates the information provided by new measurements, so that the prediction error is statistically minimized.
In the Kalman filter, three models are defined, the set of which is usually called the process model, and two phases that constitute the Kalman filtering itself: System model. Describes the evolution over time of the quantity to be estimated, expressed by means of a state vector x n+1 . The transition between states x n and x n+1 is characterized by the transition matrix A n+1 and the addition of a Gaussian white noise w n+1 with covariance matrix Q n+1 .
Measurement model. Relates the measurement vector z n to the state of the system x n through the measurement matrix H n and the addition of a Gaussian white noise v n with covariance matrix R n .
Prior model. Describes prior knowledge about the state vector at the initial time x 0 n in terms of the expected value and the covariance matrix P 0 n . Process and sensor noises are assumed to be uncorrelated.
The new value of the quantity to be estimated is predicted using the system model. For this, the estimate of the previous statex n and its covariance matrix P n are extrapolated to form the predicted state vectorx − n+1 and its covariance matrix P − n+1 , wherê Update phase. In this phase, the new state vectorx n+1 and its covariance matrix P n+1 are calculated. For this, the predicted covariance is used to calculate the Kalman gain K n+1 . The new state vectorx n+1 is calculated adding to the predicted state vectorx Water 2022, 14, 522 4 of 16 The Kalman gain is constructed to obtain minimum variance estimates. After each update, the process is repeated, taking as a starting point the new estimates of the state and of the matrix error covariance.

Cooper-Jacob Solution for the Interpretation of Aquifer Tests
The Kalman filter's implementation for the interpretation of aquifer tests was proposed based on the Cooper-Jacob solution [7] derived from the Theis Equation where: s = drawdown in the observation well located at a distance r from the pumping well [L].
where: The well function is defined as follows: This integral can be expressed as the series: Water 2022, 14, 522

of 16
Cooper and Jacob [7] noted that for small values of r and large values of t, u is small, so that in these cases drawdown can be calculated using only the first two terms of Equation (12). Therefore, Equation (9) can be approximated as: Re-writing the first term inside the parenthesis as logarithm and using the logarithm rules: which leads to the form of the Cooper-Jacob solution used in this paper for the Kalman filter implementation: In the Cooper-Jacob straight line method for the interpretation of aquifer tests, Equation (14) is converted to base 10 logs: Equation (15) plots as a straight line on semilogarithmic paper. The first step of the interpretation method consists of adjusting a straight line to the plot of drawdown (in m) versus time (in minutes) in semi-logarithmic paper; this line is projected to intersect the x-axis (s = 0).
As it can be seen in Equation (15), the slope of the adjusted straight line for a log cycle is 2.3Q 4πT . This slope is calculated as ∆s log [10] = ∆s (the value of the drawdown difference per log cycle), therefore: On the other hand, the intercepted positive value of time in the x-axis (s = 0) is designated t 0 . Substituting s = 0 and t = t 0 into Equation (15), the following is obtained: Since 2.3Q 4πT cannot be zero, for the above to hold true, it is necessary that log[ 2.25Tt 0 r 2 S ] =0, then 2.25Tt 0 r 2 S = 1. Rearranging Equation (16), the following is obtained: Thus, the procedure uses Equation (16) to solve for T and then Equation (18) to solve for S. Care must be taken to maintain consistency with the units used for the different values substituted in Equations (16) and (18). In order to avoid large errors, the critical value of u required to achieve reasonable accuracy with the Cooper and Jacob approximation is u ≤ 0.05 [4,5].

Development of Cooper-Jacob Solution for the Kalman Filter
The implementation of the Kalman filter required the mathematical work to express the Cooper-Jacob solution in the manner that this estimation method demands. Drawdown (s) and drawdown rate ( ∂s ∂t ) were selected as the state variables that compose the estimation vector, which is updated each time a new drawdown measured during the aquifer test is filtered. Therefore, it was necessary to develop the equations for the components of the transition matrix that describes the evolution of these variables through time.
First, Equation (14) is differentiated with respect to t: Substituting t 1 (time elapsed from the beginning of the aquifer test to the moment of the first drawdown measurement) and s 1 (the first drawdown measurement taken at time t 1 ) into Equation (19), the instantaneous drawdown rate at time t 1 is obtained: For B = Q 4πT , then: In an analogous manner, for t 2 (time elapsed from the beginning of the aquifer test to the moment of the second drawdown measurement) and s 2 (the second drawdown measurement taken at time t 2 ): Substituting (21) in Equation (22): Therefore, the drawdown rate at time t n+1 can be computed with the previous drawdown rate at time t n as follows: Equation (24) is valid only for t n > 0. Solving Equation (19) for s, using as integration limits s 1 , s 2 , t 1 and t 2 .
Rearranging Equation (25) to isolate s 2 and substituting (21), the following obtained: Therefore, the drawdown at time t n+1 can be computed with the previous drawdown at time t n as follows: Equation (27) is valid only for t n > 0.

Adaptation to the Kalman Filter
The propagation phase initiates with the use of Equation (4) constructed in the following manner: For Equation (5), the covariance matrix P n includes the estimate error variances of s n and ∂s n ∂t , as well as the covariance among these variables for time t n (in the first iteration it is used the prior covariance matrix). The covariance matrix Q n+1 comprises the variances for the model errors (s mod * n+1 and ∂s mod * n+1 ∂t ) and the covariance of these errors corresponding to time t n+1 . Therefore, Equation (5) is applied as: In the update phase, Equation (6) (the Kalman gain) requires the predicted covariance matrix (Equation (29)), the measurement matrix H n+1 , which has the following form: and the covariance matrix R n+1 , which includes the variance of the drawdown measurement error (s meas * n+1 ) for time t n+1 . Therefore, the Kalman gain is obtained with: Equation (7) uses the predicted state vector (Equation (28)), the measurement matrix (Equation (30)), the Kalman gain (Equation (31)) and the measurement vector. The latter is as follows: The new state vector is obtained with: The covariance matrix for the new state vector is obtained with Equation (8).
With this implementation, the model and measurement errors are provided. Therefore, it is possible to include the measurement uncertainty (associated with the technique, utilized instrument, or human error when measuring drawdown and its evolution) and evaluate how far from the model (in this case, the Cooper-Jacob solution) the real aquifer conditions are in the vicinity of the pumping well.

Optimization
To select the pair of T and S values that calibrate the Cooper-Jacob solution within the proposed procedure, we used the generalized reduced gradient (GRG), which is an algorithm for solving nonlinear optimization problems [32]. The objective is to minimize the function: whereŝ n and s C−J n are the drawdown estimates at time t n using the Kalman filter and the Cooper-Jacob solution, respectively.
The optimization procedure is subject to the following restrictions: 0.00001 ≤ S ≤ 0.001 (which corresponds to confined aquifers), T > 0, and s C−J 1 > 0. It is expected that filtered measurements will approximate the flow model estimates if measurement and model errors are adequately selected. The entire procedure was implemented in Microsoft Excel Version 2201 [33], using the solver tool for the optimization with a convergence of 0.0001 (solver The full procedure is explained in the flowchart of Figure 1. Water 2022, 13, x FOR PEER REVIEW 8 of 17 model errors are adequately selected. The entire procedure was implemented in Microsoft Excel Version 2201 [33], using the solver tool for the optimization with a convergence of 0.0001 (solver stops if the values of the relative change in the objective function are smaller than this number for the last five iterations). The full procedure is explained in the flowchart of Figure 1.

Case Study
The proposed approach was applied to confined-aquifer test data from two previously published cases. Theoretically, the assumptions employed in the development of the Cooper-Jacob solution are very approximate to real conditions for both cases.
This aquifer test was conducted in the polder "Oude Korendijk", south of Rotterdam, The Netherlands. The impermeable confining layer is composed of clay, peat, and clayey fine sand for the first 18 m below the surface. The aquifer (between 18 and 25 m below the

Case Study
The proposed approach was applied to confined-aquifer test data from two previously published cases. Theoretically, the assumptions employed in the development of the Cooper-Jacob solution are very approximate to real conditions for both cases.
This aquifer test was conducted in the polder "Oude Korendijk", south of Rotterdam, The Netherlands. The impermeable confining layer is composed of clay, peat, and clayey fine sand for the first 18 m below the surface. The aquifer (between 18 and 25 m below the surface) consists of coarse sand with some gravel. The base of the aquifer (considered impermeable) is formed by fine sandy and clayey sediments. The well screen of the pumping well was installed over the whole thickness of the aquifer, and one piezometer with a depth of 20 m was placed at 30 m from it. Since other piezometers with a depth of 30 m showed a drawdown during pumping, it could be concluded that the clay layer between 25 and 27 m is not completely impermeable. For our purposes, however, we shall assume that all the water was derived from the aquifer between 18 and 25 m, and that the base is impermeable [34].
During the execution of the aquifer test, the data in Table 1 were collected in the piezometer located 30 m from the pumping well, and the well was pumped at a constant discharge of 9.12 L 1 s −1 (or 788 m 3 d −1 ) for nearly 14 h. This case study consists of an aquifer test in a well penetrating a confined aquifer that was pumped at a uniform rate of 2500 m 3 d −1 . Time-drawdown data measured during the pumping period in an observation well 60 m away are presented in Table 2. In this case, there was no additional information about the aquifer conditions, and the constructive characteristics of the pumping and the observation well.

Results
For the implementation of the proposed procedure in the case studies, the following values were used:

The Aquifer Test "Oude Korendijk"
To construct the prior estimate vector required in Equation (28), it is necessary to select initial estimates for the drawdown and the drawdown rates. The drawdown measured at time t 1 (s meas 1 = 0.040 m) was taken as the initial estimate for the drawdown (ŝ 1 ). The initial estimate for the drawdown rate was calculated with Equation (19), resulting in: In this manner, the constructed prior estimate vector was [ 0.040 9029. 81 ].
Once we had applied the GRG nonlinear optimization for the proposed procedure, T = 510.59 m 2 d −1 and S = 0.000089 m 3 m −3 were determined.
In Figure 2, a mismatch between the simulated drawdowns and the Cooper-Jacob solution before and after the 600 min point can be observed, although both follow the general evolution of the measured data.  When using a VDME = 1.00 m 2 , T = 505.76 m 2 d −1 and S = 0.000088 m 3 m −3 are obtained. The magnitude of this considered measurement error helps to filter the measurements in a manner that allows us to obtain a better adjustment between the simulated values and the Cooper-Jacob solution (Figure 3). However, simulated drawdowns and the Cooper-Jacob solution deviate significantly from the measurements for times greater than 500 min.  For comparison purposes, in Table 3 are also presented the values for and reported in [27] obtained with traditional methods (Theis and Cooper-Jacob procedures) Figure 3. Comparison between the measured drawdowns, simulated with the Kalman filter (VDME = 1.00 m 2 ), and the Cooper-Jacob solution (the aquifer test "Oude Korendijk").
For comparison purposes, in Table 3 are also presented the values for T and S reported in [27] obtained with traditional methods (Theis and Cooper-Jacob procedures) and using an adaptative pumping test analysis based on the Kalman filter (theŞen procedure); the results of using the software AquiferWin32 version 4.00 [16] are also included. Table 3. Hydraulic parameters determined for data of the aquifer test "Oude Korendijk" using different interpretation procedures.

Parameter
Theis Procedure * In Table 3, it can be observed that the values reported in [27] for T and S differ considerably from those found with the procedure proposed in this paper. When the former parameters are used to obtain drawdowns with the Cooper-Jacob solution, significant differences are obtained with respect to the measurements, especially for larger times. As an example, in Figure 4, the Cooper-Jacob solution for T = 420 m 2 d −1 and S = 0.00016 m 3 m −3 is compared with the drawdown data of the aquifer test "Oude Ko-rendijk". A poor match is seen since the procedure for the selection of these hydraulic parameters assumed error-free measurements.
In Table 3, it can be observed that the values reported in [27] for and differ considerably from those found with the procedure proposed in this paper. When the former parameters are used to obtain drawdowns with the Cooper-Jacob solution, significant differences are obtained with respect to the measurements, especially for larger times. As an example, in Figure 4, the Cooper-Jacob solution for 420 m d and 0.00016 m m is compared with the drawdown data of the aquifer test "Oude Korendijk". A poor match is seen since the procedure for the selection of these hydraulic parameters assumed error-free measurements. The plotting of measured drawdown data showed, from the initial moments of the aquifer test, significant deviations with respect to the estimates with the Cooper-Jacob solution using the and values obtained following the traditional and the adaptative The plotting of measured drawdown data showed, from the initial moments of the aquifer test, significant deviations with respect to the estimates with the Cooper-Jacob solution using the T and S values obtained following the traditional and the adaptative methods. This result implies that these values are not representative of the actual aquifer conditions. When using the hydraulic parameters selected with the Kalman filter-based procedure, a best match between measured drawdowns and the Cooper-Jacob solution was obtained. Furthermore, increasing the value for the measurement errors produced a better adjustment between the Kalman filter estimates and the Cooper-Jacob solution; however, from 500 min after starting the aquifer test onwards, the measured drawdowns are consistently smaller than the estimates. These differences could reflect systematic errors for measurements, the non-compliance of the model assumptions, or a combination of both. Apparently, some considerable measurement errors occurred during the execution of the aquifer test, since the Kalman filter estimates and the Cooper-Jacob solution are more similar when the measurement uncertainties are increased in the proposed procedure. On the other hand, the smaller values of measured drawdowns for times greater than 500 min could reflect that some water is entering into the aquifer due to the permeability of the supposedly "confining" layers, as is stated in the description of this case study; in this case, one fundamental assumption of the selected model is not fulfilled.

The Aquifer Test of Todd and Mays
The prior estimate vector for this case study was constructed taking the drawdown measured at time t 1 (s meas 1 = 0.20 m) as the initial estimate for the drawdown (ŝ 1 ). The initial estimate for the drawdown rate was calculated as: The prior estimate vector was then [ 0. 20 2864.79 ].
Following the proposed procedure, T = 1180.43 m 2 d −1 and S = 0.00017 m 3 m −3 were determined.
The Kalman filter estimates using these values for the hydraulic parameters are very similar to the measurements and the Cooper-Jacob solution ( Figure 5), which indicates that the measurement and model errors are low for this aquifer test. The measurement taken at 60 min presents the largest error, but this is filtered during the simulation.

Discussion
The values of and obtained with the Kalman filter-based procedure were compared to those determined by other authors following previously developed approaches; the largest differences between these values were obtained for the case study with larger measurement and model errors. Following the proposed methodology for the interpretation of aquifer tests (filtering the model and measurement errors), the largest values and the smallest values were found for both case studies with respect to the methodologies that neglect these errors. According to the new values determined for these hydraulic parameters, water would flow easier, and the same amount of abstracted water would increase the drawdowns within the aquifer. Differences between the values of the hydrau- The hydraulic parameters for the different interpretations reported in Table 4 are by far more similar when compared to those corresponding to the "Oude Korendijk" case; this lower variability is related to the smaller measurement errors during the execution of the aquifer test and the high compliance of the model assumptions. For this reason, it was considered unnecessary to present the comparison between the measurements and the Cooper-Jacob solution for another pair of hydraulic parameters.

Discussion
The values of T and S obtained with the Kalman filter-based procedure were compared to those determined by other authors following previously developed approaches; the largest differences between these values were obtained for the case study with larger measurement and model errors. Following the proposed methodology for the interpretation of aquifer tests (filtering the model and measurement errors), the largest T values and the smallest S values were found for both case studies with respect to the methodologies that neglect these errors. According to the new values determined for these hydraulic parameters, water would flow easier, and the same amount of abstracted water would increase the drawdowns within the aquifer. Differences between the values of the hydraulic parameters for the different methods increase when measurement and model errors do the same.
For the aquifer test "Oude Korendijk", the largest value determined for T (510.59 m 2 d −1 ) is 49.29% larger than the smallest (342 m 2 d −1 ), while the largest value for S (0.00022 m 3 m −3 ) is 150% larger than the smallest (0.000088 m 3 m −3 ). These results reflect significant measurement and model errors for the aquifer test; it is clear that even using an extremely large value of 1.00 m 2 for the variance of the drawdown measurement error, the filtered data do not match the measured drawdowns, which reflects the considerable deviations in the aquifer conditions with respect to the proposed model.
For the case of the aquifer test of Todd and Mays, the largest value for T (1180 m 2 d −1 ) is 6.34% larger than the smallest (1110 m 2 d −1 ), while the largest value for S (0.000206 m 3 m −3 ) is 21.17% larger than the smallest (0.00017 m 3 m −3 ). In this case, measurement and model errors can be considered to be low, since a value of 0.01 m 2 for the variance of the drawdown measurement error and the proposed values for the model errors covariance matrix helps to adequately filter the measured drawdowns.
These results show the relevance of the proposed procedure for the interpretation of aquifer tests, since it provides additional information about the compliance of the model assumptions and the quality of measured drawdown data.
With the proposed methodology, drawdowns and drawdown rates can be estimated at any time of the aquifer test with their respective estimate error variances, which could be useful in cases where some aquifer test data need to be determined, are missing or seem to be suspicious.
Future research could evaluate the confidence in the execution of an aquifer test by exploring optimization methods to quantify measurement and model errors. Future developments should also be made to obtain the estimate error variances of the determined hydraulic parameters.

Conclusions
The Kalman filter-based procedure proposed in this paper represents a new alternative for the interpretation of aquifer tests in confined aquifers considering model and measurement errors. It is based on the Cooper-Jacob solution due to its relatively simple implementation in the Kalman filter estimation method.
The existing methods for the interpretation of aquifer tests produce virtually the same hydraulic parameters when measurement and model errors are low; however, the differences between the determined values of T and S increase for large measurement errors, the non-compliance of the model assumptions, or a combination of both. Neglecting these errors could lead to determining values for the hydraulic parameters that do not adequately characterize the actual aquifer conditions.
One important advantage of the proposed Kalman filter-based procedure is that it can help to quantify the measurement errors during the execution of aquifer tests and/or identify significant deviations in the aquifer conditions with respect to the model assumptions; to complement this analysis, it is fundamental to provide a description of the execution of the aquifer test (including the techniques and instruments used to measure the flow rate and drawdowns), the hydrogeological framework including the hydrostratigraphic units and the constructive characteristics of the pumping and observation wells.
Furthermore, using the hydraulic parameters selected with this procedure produces a best match between the measured drawdowns and the Cooper-Jacob solution; this result shows that the determined T and S values calibrate more accurately the selected flow model.