Influence of the spatial configuration of available data on hydraulic conductivity estimates for a geostatistical-Kalman filter method Open Access

In practice, infiltration and aquifer tests are carried out to determine the hydraulic conductivity values that characterize some areas of an aquifer. However, these field tests usually cover only a tiny portion of the entire aquifer extent. Geostatistical techniques have often been used to overcome this limitation, obtaining minimum variance estimates at unsampled sites. These estimates represent an initial spatial configuration of the hydraulic parameter for the construction of conceptual models, the execution of water balances, or as reference fields for groundwater model calibration processes.

Geostatistical methods have been compared with different interpolation methods. For example, Panagopoulos et al. (2006) compared ordinary kriging, inverse distance, and Thiessen polygons; the geostatistical kriging method gave the best results. Similarly, Ghosh & Pekkat (2019) compared kriging, inverse distance weighted, natural neighbor, spline, and trend interpolation methods for the spatial prediction of hydraulic conductivity; they got the best results with kriging even for cases with a highly heterogeneous medium with the additional benefit of providing a measure of the estimate uncertainty, an important feature of geostatistical techniques with respect to traditional interpolation methods. For these cases, a good agreement between the obtained errors and the estimate error variances gives reliability to the selected variogram model.

The hydraulic conductivity can be estimated within an aquifer using data obtained from field or laboratory tests (Ahmad et al. 2021; Chandel et al. 2022). However, these data are usually scarce and characterize only small portions of the total area of an aquifer; additionally, hydraulic conductivity varies by several orders of magnitude in small regions. For these reasons, the spatial configurations of K obtained with traditional interpolation methods present a high uncertainty in most cases. In this sense, bivariate geostatistical techniques (using a more sampled secondary variable) have represented an alternative to improve hydraulic conductivity (K) or transmissivity (T) estimates reducing the estimation uncertainties.

Freixas et al. (2017) revealed preferential flow paths by analyzing cokriging transmissivity estimates using the storage coefficient as a covariate. Also, Patriarche et al. (2005) used two approaches to estimate hydraulic conductivity. In the first, cokriging was used to estimate K (direct approach), and in the second, cokriging was used to estimate transmissivity and K was retrieved dividing T by the well screen length (indirect approach); for both, the specific capacity was the covariate. The best results were obtained by the indirect approach: log (⁠⁠) data, and dense T sampling. Ahmed & de Marsily (1987) also used specific capacity as a secondary variable. Frohlich & Kelly (1985), Kupfersberger et al. (1992), El Idrysy & De Smedt (2007), Xu et al. (2017), and Chen et al. (2018) included electrical resistivity as secondary data to get K or T spatial estimations when electrical resistivity and the hydraulic parameter have a good correlation, and the electrical resistivity data is densely sampled. Other secondary variables have been used to estimate T using geostatistical methods, for example, seismic data (Rubin et al. 1992), and slope of the water table (El Idrysy & de Smedt 2007).

As far as we know, Ahmed & de Marsily (1989) were the first to test the use of transmissivity and piezometric measurements simultaneously to estimate transmissivity (T) with cokriging and compared the T estimates with those of a kriging approach. The hydraulic head (HH) and the transmissivity natural logarithm (lnT) were assumed to be related through an equation. Thus, the covariances of lnT, HH, and their cross-covariance were related through a differential equation. A set of Gaussian covariance models that satisfied these differential equations was proposed. The authors tested their approach using a hypothetical example based on 2,000 km² of a chalk aquifer in northern France. A conditional geostatistical simulation of T was generated based on 72 actual transmissivity measurements. Thirty-five synthetic T and 144 HH measurements were generated using a flow model with the simulated T. The model did not include sources or sinks. For the geostatistical analysis, the HH trend was removed using a simple flow model; the fitted ln T covariance was isotropic, and the covariance of the HH residuals and the cross-covariance of HH and T were anisotropic. For the kriging T estimates, simple kriging was used, whilst residual cokriging was used for the bivariate estimates. The results indicate that the errors of the cokriging estimates were smaller than the kriging ones. However, the mean squared error (MSE) of cokriging was larger than that of kriging (5.79 and 2.4 m²/d, respectively), indicating that the cokriging variance underestimated the errors.

Most work undertaken to estimate K or T using piezometric data has been done to solve the inverse problem for groundwater flow and transport models (Kitanidis & Vomvoris 1983; Xu et al. 2013; Erdal & Cirpka 2016) maybe because the models establish a natural relation between K or T and HH. In the context of stochastic groundwater modeling, ensemble Kalman filter methods have been used to estimate K using HH or drawdown measurements (Panzeri et al. 2015; Xu & Gómez-Hernández 2015; Carrera et al. 2020). Other approaches to estimate saturated hydraulic conductivity in large areas of an aquifer consist of using artificial intelligence algorithms (Gorgij et al. 2018; Yusefzadeh & Nadiri 2021; Singh et al. 2022). However, to our knowledge, none of them perform a sensitivity analysis of the number and spatial distribution of the available data.

Júnez-Ferreira et al. (2020) took up the use of bivariate geostatistical-based techniques to estimate hydraulic parameters within an aquifer. They present univariate, bivariate, and bivariate space-time geostatistical-Kalman filter-based approaches (similar to kriging, cokriging and spatio-temporal kriging) to estimate K in heterogeneous aquifers. For the bivariate case, data were used. The correlation structures between data of the principal (⁠⁠) and the auxiliary (⁠⁠) variable were determined through geostatistical analyses and the estimates were obtained through a static Kalman filter (that does not use the model equations). Data from a synthetic numerical flow model were used to test the method. A fixed number of spatially uniformly distributed K and data were used to estimate K with the static Kalman filter.

Several authors have pointed out the importance of using reliable variogram estimates to reduce the estimation error in kriging interpolations (for example, Gumiere et al. 2014; Oliver & Webster 2014). The number and spatial locations of available data directly influence the selected variogram model and its parameters and, therefore, the variable estimate. But, small changes in the model function and its parameters can significantly affect the kriging error variances, which are essential to predict the estimation error magnitude. Oliver & Webster (2014) recommend using 100–150 data to obtain reliable variogram estimates. However, this recommended highly dense sampling is not a common situation in many real cases, specifically for hydraulic conductivity. Also, as far as we know, general criteria on the number of data needed to get reliable estimates using cokriging are not available. Therefore, evaluating the estimation method's performance under different numbers and spatial configurations of the available data is advisable. Even though geostatistical methods can give reasonable K estimates, it is essential to evaluate the model performance when the estimates are used for numerical modeling purposes.

This paper presents a sensitivity analysis for the univariate and the bivariate estimation methods proposed by Júnez-Ferreira et al. (2020), evaluating the hydraulic conductivity estimates for different numbers and spatial locations of the available data for the principal and auxiliary variables. These estimates are compared with the known hydraulic conductivity values of a synthetic numerical flow model. This evaluation aims to support the design of optimal sampling strategies (for the measurement of hydraulic head and the execution of pumping tests) oriented to obtain reliable hydraulic conductivity estimates. These actions, among other benefits, could reduce calibration efforts by providing adequate reference fields and increasing the knowledge of the correlation structure for field variables.

The remaining sections of this paper are structured as follows. In the Materials and Methods section, the geostatistical theory and static Kalman filter used in this paper are presented, together with the structure of the covariance matrices for the univariate and bivariate cases. The Case Study section includes the characteristics of the synthetic numerical flow model employed and the datasets used for the estimation of K and its comparisons; it also includes the error statistics used for these evaluations. In the Results and Discussion section, the parameters of the selected variogram models, the evaluation of the K estimates, and simulations of the numerical flow model for the univariate and bivariate cases are shown and discussed. Finally, Conclusions are presented.

The influence of the data number and spatial locations on the hydraulic conductivity estimation was evaluated for univariate and bivariate geostatistical–Kalman filter estimation approaches. The evaluation was done using data generated from a previously built transient numerical flow model. Different number, spatial arrays, and combinations of K and simulated were used for estimating K. Different variogram models for the univariate and bivariate cases are obtained from geostatistical analyses. They are used to obtain the corresponding covariance matrices that the static Kalman filter employs for estimating K. The hydraulic conductivity was estimated by two different approaches:

• - Univariate estimation (based on the spatial distribution of hydraulic conductivity data only);

• - Bivariate or cross estimation (hydraulic conductivity as the primary variable and hydraulic head for the last simulation as the secondary variable).

A description of the geostatistical procedure to derive the univariate, and bivariate covariance matrices is presented in the following subsections. Relevant aspects of the static Kalman filter theory and its implementation for the proposed methodology are also included.

Geostatistical analyses were done to get the spatial and cross structures of covariance matrices derived from variograms.

For the recursive implementation of these equations, a prior state estimate and a prior covariance matrix estimate are required; the latter is obtained from a geostatistical analysis for the case study. After each update, the process is repeated, starting from the last state estimate and its error covariance matrix (Briseño-Ruiz 2012).

For each new sampling data, the Kalman filter modifies the state vector estimate and the covariance matrix estimate using equations 6, 7, and 8. For the univariate case, the available hydraulic conductivity data are used to update the estimate, and for the bivariate case, the hydraulic conductivity and hydraulic head data are used.

In this paper, the state vector is ⁠, for the univariate case, and ⁠, for the bivariate case, where has components ⁠, is a given position, and Y is hydraulic conductivity or its natural logarithm ⁠, depending on the test case in question (they are explained in the Case Study section below); HH has components . Y length is ⁠, and length is ⁠, so the total length is ⁠, which is also the number of files and rows of the covariance matrix ⁠. Finally, one measurement is processed at a time, so is a zero-row vector but at the measurement position at which its value is 1. For the bivariate case, even when the vector includes and ⁠, 's estimate is the only one used. The corresponding subvector of the prior state vector equals the available data mean value of Y or ⁠.

The different error statistics used to evaluate the performance of each analyzed case are given below:

Different sets of K and data (corresponding to the final simulation time results) were used to estimate K and in the modeling grid nodes. The resulting K and estimates were compared with the reference K and values (Figure 1(a) and 1(b), respectively). A flow direction from South to North can be seen imposed by the defined boundary conditions, with lower hydraulic gradients in the North where the larger K values are found. Furthermore, equipotential lines are closer to straight lines in the North half while a zig-zag type behavior is produced in the contrast of low K values in the South.

For the univariate case, subsets of the 400 hydraulic conductivity reference values were used as ‘data’. The analyzed cases use 20, 40, and 60 K data, selected arbitrarily but seeking a spatial distribution close to homogeneous through the entire aquifer, and 20 K data distributed: (a) in the Northern half of the aquifer, and (b) in the Western half of the aquifer. Additionally, six analogous cases using natural logarithm of K data were analyzed.

For the bivariate case, the hydraulic conductivity was estimated using 20, 40, and 60 K data, or their natural logarithm, distributed through the entire aquifer, in combination with 60 and 400 data for the last simulation time (used as the secondary variable), and 20 K data, or their natural logarithm, distributed in the Northern and Western halves, combined with 60 or 400 data. The 400 data cover the entire estimation grid, and the 60 data locations are the same as for the 60 K data.

 The state vector length is N = 400, and the initial state is the average value of the available data (⁠ or ⁠) for each case. Then, the covariance matrix for each case has 400 rows (structured as shown in Equation (9)), its components are calculated using the selected variogram model for the ‘available’ hydraulic conductivity data for each case. These variograms are shown in Table 1; note that larger sill values are registered when K values are distributed through the entire aquifer extent.

Estimates on the finite-differences grid nodes for all cases were obtained using the Kalman filter and compared with the reference values (the 400 K values of the numerical flow model), and the errors were calculated as the difference between them. Besides the minimum, maximum, and mean error, the RMSE, SMSE, and Mean Absolute Percentage Error (MAPE) are reported in Table 2. The best results for each sampling case are highlighted in bold (in Table 2), with regard to the lowest errors and the agreement between these errors and the estimate error variances. For the cases of 20, 40, and 60 K data distributed through the entire aquifer, the best SMSE values (closer to one) correspond to the lowest variances which gives more confidence to those estimates. For the samples distributed at north and west, the best SMSE values correspond to the largest variances. Although the lower MAPE value is obtained for the log-transformed 60 K data option, the overestimation of the estimate error variances is reflected in the corresponding low SMSE value which reduces the reliability of the selected variogram model.

The best estimates were obtained, as expected, when 60 K data distributed through the entire aquifer extent were available, in particular, with K data without transformation. The RMSE value is the best since it is closer to zero than for the other cases; even though the value for SMSE reflects an overestimation.

For the 20 K data estimates, K is overestimated in the North Central and Southern half regions (the lowest K reference values are located in the latter). The higher estimated error variances are in the Southeast and South-Central regions with values above 4 m²/d².

For the 40 K data case, K is overestimated in the northern half and the middle regions; the South is better characterized than the 20 K data estimates. Low estimate error variances (between 1.01 and 2.00 m²/d²) are present in a more extended area; on the other hand, in the Southeast region, high values, above 4 m²/d², are present.

For the 20 K (N) data case, K is overestimated in the North Central and in the Southern half, where all the estimates fall between 6.01 and 7.00 m/d due to the lack of data in these regions. The estimated error variances in the Northern half are similar to those of the 40 K data case. However, the Southern half presents a highly extended area with values between 4.01 and 5.00 m²/d².

The different hydraulic conductivity estimates were used as input for the K values in the numerical flow model and the simulated results were compared with those of the original model. These ‘errors’ were used to evaluate how well the hydraulic conductivity estimates allow the flow mechanisms of the original model to be reproduced (Table 3; the best simulations for each sampling case are highlighted in bold). The simulations with the lowest errors correspond to the estimation using 60 K data, differences with simulations of the original model (400 spatially distributed data for 6 time-steps) were minimal. As long as fewer K data were used for the estimation, larger errors were obtained. The maximum errors correspond to the case of K data grouped in the North followed by the case of K data grouped in the West.

The MAPE shows that for all cases the model error represents a minimal fraction of the data at each node; this means that all the estimations are robust enough to be used in the numerical flow model which has low sensitivity.

The state vector length for the bivariate case is twice the size of that for the univariate case because it now includes the 400 s of the final simulation time. The prior estimate values for each subvector are the average data values ‘available’ for K or ⁠, respectively. The K and covariance matrices and the cross-covariance matrix are derived using the selected variogram models for each case's K and ‘available’ data (Tables 4 and 5). Each matrix is constructed of 400 rows and 400 columns (as stated in Equations (9, 10 and 11) for the bivariate case).

Cross-validation was used to evaluate the best variogram model fit from all the analyzed models. As can be seen in Table 6, when the K data are distributed through the entire aquifer extent, the increase in the number of data is almost unnoticeable in the value of the variogram parameters except for the range; on the other hand, changes in the number of K data produce slightly more significant changes in the nugget and sill values. The impact of the number of data in the variogram parameters is larger when the K data are grouped in the North and the West.

Tables 4 and 5 show the correlation between the or and data. The higher correlations on average correspond to the cases when 20 are grouped in the North, followed closely by the cases where 20 K data are grouped in the West. The lower average correlation was for the cases where 40 K data are distributed throughout the extent of the aquifer.

The estimate errors for the bivariate cases are shown in Table 6; the best options for the different number of available K data are highlighted in bold. For the 20 K data distributed through the entire aquifer extent, the lower estimate errors were obtained for one of the cases with higher correlations between K and data (the 20 – 60 data case). The second lowest variances for these cases were obtained for the case of 20 K and 60 HH data, with a better agreement between the calculated errors and the estimate error variances, which makes these estimations more reliable. Using the data over the 20 K univariate cases reduces the RMSE for K by up to 16.96%.

For the bivariate cases that use 40 K data, the best estimates were obtained for 40 K – 60 data followed closely by the 40 K – 400 data (both with lower correlations among the 40 K cases). A marginal improvement in estimations is noticed compared to the univariate cases (the RMSE for K is reduced by up to 4.71% for the lnK cases). In both bivariate cases (without transformation) the SMSE values are also very similar, however the 40 K – 60 HH data are preferable since less data are required for similar estimations.

When 60 K data are available, the best estimates were obtained for the univariate case without transformation; however, the 60 K – 60 and the 60K– 400 HH data cases are similar with even better SMSE values. The 60 K – 60 HH data option was selected as the best since a lower number of HH data gives practically the same results. It is important to note that the RMSE for the 60 lnK univariate case is 1.28% lower than that obtained for the best bivariate case.

For the bivariate cases of 20, 40 and 60 K data distributed through the aquifer, the RMSE values for the best options are similar to the respective univariate cases, however the SMSE values for the former cases are closer to one; additionally the corresponding estimate error variances are lower which provides more confidence to the estimates obtained for the bivariate cases. The MAPE value is improved substantially for the 20 and 40 K data. The contribution of the secondary variable is reduced for a larger amount of data for the primary variable.

With respect to the K data grouped at in the North or West, larger variances are obtained for the bivariate cases, but according to the SMSE there is a better correspondence with the estimate errors. Furthermore, the calculated errors are lower for the bivariate cases with a notable reduction of the MAPE for the best bivariate case of K data distributed in the north, and for the best bivariate case of K data distributed in the west. These results clearly show the major contribution of the secondary variable in those zones where the primary variable has been poorly sampled.

In general, the best estimates were obtained for the bivariate cases with higher correlations, for the 20 K data cases. The worst estimates were produced for the bivariate cases of 20 K data grouped in the North (N) and West (W), with similar correlations for both (Table 5). The best estimates for these cases are obtained when using 20 (N) – 400 and the 20 (W) – 400 data (with the second lower correlations among their respective cases). With respect to the univariate cases using 20 K data grouped in the North (N) and West (W), the RMSE for K is reduced for their respective bivariate cases.

For the 20 -60 data bivariate estimates of K, this parameter is better estimated in the Northwest and Southern regions when compared to the univariate case. Furthermore, the highest estimate error variances (values above 4 m²/d²) can be found in more extended areas in the Southeast and South-Central regions. In general, the estimate error variances are higher for the whole aquifer.

The bivariate estimates of K for the 40 -60 data are very similar to the 40 K univariate estimates, with an overestimation in the northern half and the middle regions. Low estimate error variances (between 1.01 and 2.00 m²/d²) are found over a more extended area when compared to the univariate case, high values above 4 m²/d² can be found in a small part of the Southeast.

Estimates of the bivariate – 400 data characterize the North Central and the middle part of the aquifer in a better manner; on the other hand, K estimates are now between 6.01 and 7.00 m/d in a more reduced area of the South. The estimated error variances are much lower in the Northern region; however, the extreme South presents an extended area with values between 6.01 and 7.00 m²/d², which are higher than the univariate case.

For the cases with 20 K data distributed within the entire aquifer extent (with or without using a secondary variable), the minimum RMSE do not correspond to the case where the SMSE values are closer to one; the variance values are large for the cases of the minimum RMSE. On the other hand, for the univariate and bivariate cases using 40 and 60 K data, the minimum RMSE correspond to the best agreements between the calculated errors and the estimate error variances obtained with the selected model (SMSE values closer to one).

For the cases of 20 K data distributed in the North or in the West, the minimum RMSE correspond to the cases of the second best SMSE values.

After the bivariate K estimates were used to simulate the hydraulic head evolution, the model hydraulic head errors were calculated (Table 7). The best simulations for each sampling case are highlighted in bold. The best results for 20 K data distributed through the entire aquifer extent correspond to the 20 – 60 and 20 – 400 cases (the error statistics are very similar) which correspond to the best K estimates and the highest correlations of these cases. The improvement in the simulations concerning the univariate case is notable (the RMSE for the hydraulic head simulation is reduced by up to 56.67%). The 20 – 60 case is considered better since it demands less data for the secondary variable.

The best simulation for the bivariate cases using 40 K data corresponds to the 40 – 60 case. Significant improvements in simulation results are obtained compared to those in the univariate case (up to 41.67% reduction of the RMSE for hydraulic head simulation), but lower than the best for the 20 K cases.

For the bivariate cases using 60 K data, the best simulations were obtained for the 60 – 60 (with the highest correlation). The advantage of the bivariate cases compared to the simulations with the univariate case K estimates is also clear (the RMSE for hydraulic head simulation is reduced by up to 53.33%).

As for the univariate cases, the MAPE shows that for all the cases the calculated error represents a minimal fraction of the data at each node, which corroborates the robustness of the estimates to be used in the numerical flow model. The major reduction in the MAPE corresponds to the simulation obtained for the best bivariate case with 20 K data grouped in the west (it is reduced 0.05% with respect to the corresponding best univariate case).

Although the log-transformed K data do not help to obtain the best K estimates for the bivariate cases, they provide the best results in the numerical flow model simulations for all the bivariate cases of K data distributed through the entire aquifer extent. This transformation favors obtaining similar absolute values of the bivariate variogram parameters for the and the sills and a larger value of the – sill component. With these similar values, a higher weight is given to the data in the estimation of K (without logarithmic transformation, the largest sill value corresponds to the univariate K sill component), allowing more advantage to be taken of the information provided by the secondary variable.

For the cases when the available K data are located in the northern portion of the model domain, the best simulations are obtained for the 20 K (N) – 60 data (which coincides with the best estimates and the second lowest correlation of these cases). The best simulations for the cases where the K data are grouped in the West portion of the domain correspond to the 20 K (W) – 60 data (it has the second-best correlation of these cases). Bivariate cases for data grouped in the North and West offer better simulations than their respective univariate cases, due to the additional information provided by the data in areas without K data. Unlike the K data distributed through the entire aquifer extent cases, the worst simulations among these bivariate cases are obtained when the logarithmic transformation is used. Using 60 instead of 400 data seems to be a better option according to simulation results. The RMSE for hydraulic head simulation is reduced by up to 14.14 and 78.21% for the bivariate K (N) and K (W) cases, respectively.

As expected, K estimates are improved as the spatial coverage of K increases, reducing the contribution of data. In the Kalman filter implementation, the use of K data rapidly reduces the estimate error variance on estimation points; therefore, for bivariate cases, information provided by data is highly redundant in the vicinity of K sampled sites. Therefore, more benefit is obtained from the information provided by the secondary variable when there are extensive areas without data for the primary variable, as occurs with the case of K data grouped in the West; since there is no data for K in the East, the larger the number of data, then the estimates are better in that zone.

However, care must be taken in obtaining the cross-variogram and selecting the initial state to maximize the advantage of using a secondary variable. In this sense, prior to implementing these bivariate techniques it is necessary to evaluate correlation between variables and planning for an adequate sampling scheme that considers hydrogeological knowledge of the study area.

Ahmed & de Marsily (1989) presented the only previous geostatistical study we know of in which cokriging is used to estimate T using HH data (not in the inverse problem context), and the results are compared with kriging T estimates. In this work, we use a static Kalman filter to estimate K, which is closely related to T, and we also use HH data. Both methods are related because they are best linear unbiased estimation methods, and, in both, the spatial structure model is obtained through a geostatistical analysis. Our results confirm that cokriging cross-validation estimate errors are smaller than the kriging errors for the corresponding case; however, in contrast to their results, our results had SMSE values closer to one than the kriging values were. On the other hand, we point out that Ahmed & de Marsily (1989) did not test the T estimate's sensitivity to the data's number and locations.

Other papers estimate hydraulic conductivity using hydraulic head data in the inverse problem context. However, we note that an inverse problem objective is different from ours, so we do not think comparing our results with a paper of that nature would provide insights into our analysis. Recent works that estimate hydraulic conductivity using artificial intelligence also do not perform a sensitivity analysis of the estimates to the number of data and their spatial distribution.

The results show that the use of HH data as a secondary variable in the evaluated geostatistical–Kalman filter estimation method helps to obtain low estimate error variances in the estimation of K; furthermore, from a geostatistical point of view, these results are more reliable since they have a better agreement with the calculated errors from the cross validation of the selected variogram model. However, it is noticed that the contribution of the HH data (as secondary variable) is reduced when the number and coverage of available K data increases. A major contribution of the HH data is achieved in those areas of the aquifer where this variable has been measured and K has been poorly sampled.

The differences in the simulated values of HH for the different spatial configurations of K obtained with the geostatistical–Kalman filter method for the case study are minimal, which implies the robustness of the estimation method to obtain reference fields useful for a groundwater model calibration process.

It is also worth mentioning that the method proposed in Júnez-Ferreira et al. (2020) is faster than the kriging and cokriging methods because the Kalman filter allows processing one data at a time, avoiding inverting a matrix. Processing the test case with 60 K and 400 HH data (the one with the most data) takes a few seconds on a notebook computer with an Intel Core(TM) i76500U CPU (2.5 GHz, 2.6 GHz) and 16GB RAM.

Univariate geostatistical techniques, such as ordinary kriging, represent a standard tool in hydrogeology to obtain low uncertainty estimates of hydraulic conductivity; however, more recently, the use of auxiliary variables has been promoted to improve the results. In this paper, the influence in the estimation of hydraulic conductivity of the number and spatial location of data was evaluated for univariate and bivariate Kalman filter–geostatistical approaches, and the errors of the hydraulic heads simulated with the resulting K estimates.

The use of statistics that measure differences between estimated and reference values is perhaps the most usual manner to evaluate estimation techniques; however, in most field situations, there are no data to validate the estimation results, making synthetic cases an excellent alternative to do the evaluation when introducing and validating new methods.

In the case study presented, it was found that bivariate estimates take advantage of the correlation between the primary and the secondary variables (⁠ and ⁠, respectively) to outperform the univariate estimates (the improvement in the estimation for bivariate cases is reduced when a better sampling scheme for the primary variable is followed). Our analysis shows that hydraulic head can provide valuable information where hydraulic conductivity was not sampled.

For those areas where only is available, the bivariate estimated K outperforms the univariate cases, but the obtained estimate error variances are larger than in the univariate case. The larger sill of the bivariate K model variograms can explain these large variances.

The variation in the number of data (the primary variable) produces more changes in the bivariate variogram than the variation of the number of data (the secondary variable). This is because ‘available’ data allow the estimate error variance for and to be reduced to close-to-zero values, and the influence of the secondary variable (⁠⁠) data on the estimates is insignificant.

For the bivariate cases with a uniform spatial distribution of the K and data, a higher correlation between K and data does not guarantee better estimates of K. When K is uniformly distributed over the estimation area, the data contributes more to improving the K estimates if the largest sill value of the bivariate variogram model corresponds to the univariate component of ⁠. Also, when K is not sampled for the entire estimation area, a higher correlation between the primary and the secondary variable and a good selection of the initial state vector are more relevant.

The best model simulation results for the bivariate cases with data distributed through the entire aquifer extent are obtained using (corresponding to the highest correlations); however, this option does not always offer the best estimates. The cases using 20 K data grouped in the Western half of the domain only (20 and 400 data) are very close to those obtained with the best simulation of the bivariate cases using 20 K data uniformly distributed (20 and 60 data). These results corroborate the advantage of using a secondary variable to improve the estimations of the primary variable at unsampled sites. For the case study, the sampling of K in a fringe from North to South (aligned with the flow direction) takes more advantage of the correlation between K and HH in comparison to the case where K is sampled in a fringe from West to East (perpendicular to the flow direction).

One of the major limitations of this study is that the proposed synthetic case is simple. It consists of a single layer isotropic medium with a constant storage coefficient, and a simple geometry. Also, the measurement errors for K and HH were neglected.

One important improvement of the estimation method would be the inclusion of additional auxiliary variables such as electrical resistivity, specific capacity and water level slopes.

For the case study, the sampling of K data oriented in the flow direction takes more advantage of the HH data located in unsampled areas of the aquifer; this result deserves a deeper analysis in order to optimize the characterization of the hydraulics parameters of an aquifer. Another interesting topic is the estimation of anisotropic K fields.

Future work may include strategies for sampling design (including the hydrogeologic framework) to improve estimates of within an aquifer, considering uncertainty in ‘available’ data, and more complex cases.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.