ABSTRACT

Groundwater contamination by nitrate within an unconfined sandy aquifer was mapped using a Bayesian Data Fusion (BDF) framework. Groundwater monitoring data was therefore combined with a statistical groundwater contamination model. In a first step, nitrate concentrations, measured at 99 monitoring stations irregularly distributed within the study area, were spatialized using ordinary kriging. Secondly, a statistical regression tree model of nitrate contamination in groundwater was constructed using land use, depth to the water table, altitude and slope as predictor variables. This allowed the construction of a regression tree based contamination map. In a third step, BDF was used to combine optimally the kriged nitrate contamination map with the regression tree based model into one single map, thereby weighing the kriged and regression tree based contamination maps in terms of their estimation uncertainty. It is shown that BDF allows integrating different sources of information about contamination in a final map, allowing quantifying the expected value and variance of the nitrate contamination estimation. It is also shown that the uncertainty in the final map is smaller than the uncertainty from the kriged or regression tree based contamination map.

1. Introduction

Assessing the quality of groundwater is a prerequisite for designing sustainable water management strategies e.g. when implementing the Water Framework Directive [1]. Groundwater quality is, however, a spatially distributed attribute. It is generally accepted that knowledge of the spatial distribution of this attribute allows the design of site specific protection and remediation measures. Therefore robust and validated techniques are needed to map the spatial distribution of groundwater quality within the groundwater body continuum. Unfortunately, properties related to groundwater quality, such as nitrate concentration, can only directly be measured at the local scale. The mapping of groundwater quality within the water body continuum will therefore often build on data interpolation or prediction of properties related to groundwater quality through deterministic or stochastic models, and the adopted technique will impact the results of the final assessment.

Point data can be spatialized with traditional interpolation tools such as Voronoi tessellation, inverse distance weighting or kriging. References [2,3] for instance used ordinary kriging for mapping the spatial variability of nitrate concentrations in a shallow water body. Reference

[4] used Gaussian simulation techniques to introduce local uncertainty for mapping nitrate concentrations within the groundwater bodies of the Po Valley (Italy). Among modern interpolation techniques, the Bayesian Maximum Entropy (BME) framework was developed recently by [5] as a generalization of classical techniques. Reference [6] used BME as a formal spatial modeling framework allowing introducing the temporal variability of groundwater contamination and the sampling rate in a spatial map.

Unfortunately, the locations where groundwater quality properties are measured are typically scarce and sparsely distributed over space. As a result, large uncertainties about the variable of interest arise in the mapping process, especially at prediction points located far away from monitoring stations. Given that the quality of the interpolation decreases rapidly with the distance to the monitoring station, the mapping of a variable of interest through classical geostatistical interpolation techniques will have limited practical utility in regions where a densely distributed monitoring network is not available. Unfortunately, this is rather the rule than the exception, and therefore alternative robust and validated methods are needed to provide more comprehensive and accurate information about the groundwater quality. Data fusion is such a method that can be used to combine optimally various sources of information about groundwater quality in a consistent and accurate model prediction.

Data fusion techniques have already been used in various application fields, including for example biometrics verification systems (e.g. [7,8]), surveillance systems (e.g. [9]), robotics (e.g. [10]), medical imagery (e.g. [11]) or military/civil engineering (e.g. [12,13]). Data fusion has also important applications in classification of remote sensing images (e.g. [14]) and in environmental modeling [15].

Among different data fusion techniques, a Bayesian Data Fusion (BDF) approach was recently proposed by [15]. It was especially designed for spatial predictions problems and provides a consistent framework of fusing an arbitrary large number of information sources that are related to a same variable of interest in order to provide a unique spatial prediction. The main advantage of a Bayesian approach is to put the problem of data fusion into a clear probabilistic framework. Since this original paper, the general method has been applied to various case studies of environmental sciences; remote sensing [16, 17], hydrology [18,19] or air pollution [20]. Similarly to these applications, the present paper relies on some specific assumptions with the aim of spatial mapping of groundwater quality. To the best of our knowledge, BDF has never been used for mapping groundwater pollution by nitrate.

In this paper, we used the BDF approach as a formal framework to map groundwater contamination by nitrates. Groundwater contamination by nitrate remains a critical water quality issue for many groundwater bodies all over the world. The approach was illustrated for mapping groundwater contamination in the water body of the Brusselian sands located in an unconfined sandy aquifer in the centre of Belgium. This water body is considerably affected by nitrate pollution.

Groundwater nitrate concentrations were measured at different monitoring stations irregularly distributed in the study area. In a first step, these point data were interpolated on a grid covering the whole study area by a classical (ordinary kriging) geostatistical technique. With this interpolation method, the uncertainty on the prediction depends partly on the respective geometry of the monitoring station and the estimation point. As a consequence, the uncertainty on the predictions is larger in parts of the study area where data are scarce, making these estimates useless at some points on the map. Secondly, nitrate concentrations were predicted by a statistical model all over the study area. With this model, the prediction uncertainty depends on the uncertainty of underlying predictor variables such as depth to the water table, land use or altitude, along with the model uncertainty. In a third step, nitrate concentrations estimated by the interpolation method were combined with those obtained from the statistical model into a single prediction. In this step, BDF is used. The BDF prediction together with its estimated uncertainty was further spatially mapped.

2. Study Area and Data

This study focuses on the unconfined sandy groundwater body located in the Brusselian aquifer in the center part of Belgium. The aquifer has a surface area of 965 km² and is of primary importance for drinking water supply. This unconfined aquifer is located in Tertiary sands and is overlaid by a Quaternary loess layer of variable thickness (0 to 15 m). The Brusselian sands outcrop mainly in the valleys where sandy and sandy loam soils develop. Transmissivity of the aquifer varies from 2.9 × 10^–5 to 1.2 × 10^{–2 m2}/s and its permeability varies from 1.4 × 10^–6 to 6 × 10^{–3 m/s [21].}

This aquifer is characterized by both the presence of intensive arable cropping and intense urban pressure. The 1:10.000 land use map with 65 land use classes of the Walloon Region was provided by the regional administration and depicts the situation of 2005 [22]. The land use is highly fragmented. Typical land uses are urban (generally located in the valleys; about 17% of land use in the study area), grassland and forests (found on valley slopes; about 13% and 10%, respectively), and arable land, mainly wheat, sugar beet, maize and barley (found on loamy soils on the plateau; about 51% of land use).

The depth to the water table was calculated by subtracting the piezometry value from a 30 m resolution digital elevation model (DEM) which was furnished by the regional administration. The piezometry map was calculated by interpolating the water table levels measured in 1984 using ordinary kriging. The calculation of depth to the water table varies from 0 to more than 45 m with a mean value of 10.5 m.

The groundwater nitrate concentration data used in this study was recently collected over 99 monitoring stations in January and February 2009. These monitoring stations are wells, galleries, drains and springs. A one-way analysis of variance was performed for comparing the nitrate concentrations measured in the different types of monitoring stations and no significant differences were detected in the mean concentrations of the different groups at a significance level of 0.05. The nitrate concentrations show a wide spatial variability in the study area, with values going from 6.9 up to 93.4 mg NO₃/L. The regional mean and standard deviation are respectively equal to 45.6 and 17.2 mg NO₃/L (Table 1). A Lilliefors test was performed using the function implemented in the Statistics Toolbox^TM of Matlab^TM. This test confirmed the Gaussian shape of the distribution of the measured nitrate concentrations at a significance level of 0.01.

3. Tools and Methods

3.1. Kriging

Kriging is a group of stochastic prediction techniques widely used in geostatistics to interpolate the value of a random field (e.g., the groundwater nitrate concentration) at an unobserved location, based on a linear combination of observed values at nearby locations. Kriging incorporates the spatial dependence of the data in its estimation process through a variogram or a covariance function. The variogram function yields the average dissimilarity between locations separated by different intervals of distances. Kriging is known to provide a linear predictor that corresponds to the Best Linear Unbiased Predictor (BLUP) in the least squares sense. Additionally, it is the best possible predictor when random field is assumed to be multivariate Gaussian. Ordinary kriging is the most commonly used type of kriging. It assumes a constant but unknown mean and enough observations to reliably estimate the variogram.

3.2. Regression Tree

Regression tree modeling is an explanatory technique relying on a process known as binary recursive partitioning. Regression trees became popular in environmental sciences in the early nineties (e.g. [23-25]). The algorithm identifies which of the variables explains most of the variance in the response variable, then determines the threshold value of the explanatory variable that best partitions the variance in the response such that it minimizes the sum of the squared deviations from the mean in the two groups. The process is repeated for each new branch until there is no residual explanatory power, according to the limitations imposed by the user. Suppose that predicttor variables x₁, x₂, ···, x_N and the response variable y are organized as columns in a table. The database table is sorted by the column of the first variable (x₁). Then the table is split into two parts (called left and right branch) until the number of samples in the branch to split is lower than a user defined threshold. For each possible split position, the two partitions are compared based on the reduction in non-homogeneity

(1)

that they provide. The non-homogeneity in a group of samples is measured by computing deviations and is defined as

(2)

where is the mean value across all observations y_i. Each partition generates left (D_L) and right (D_R) deviance values. This process is repeated for each of the N predicttor variables. The partition that maximizes the change in deviance ΔD is the partition to choose. Each of the branches obtained after partitioning is partitioned again using the same method.

Regression trees are often used to see whether complex interactions between explanatory variables exists and to identify which one of the predictors have the most important effect on the dependent variable [26]. Major advantages of regression tree models are that 1) they are nonparametric and, hence, Gaussian distribution assumption of predictor variables does not need to be satisfied, 2) they can incorporate categorical data, and 3) they allow possibly complex interactions between the predictor variables to be represented without assumptions of linearity. Furthermore, while multiple linear regression identifies global relationships in the data set, regression tree are able to identify local relationships [27].

In this study, a regression tree was developed for predicting nitrate concentrations as a function of the predicttors listed in Table 2. The regression tree model was

applied to the dataset to predict nitrate concentrations all over the study area at the nodes of a regular 50 m × 50 m grid.

Regression tree algorithms are implemented in various software packages such as Matlab^TM, SAS^TM, R. The “classregtree” function implemented in the Matlab Statistical Toolbox^TM was used in this study.

3.3. Bayesian Data Fusion (BDF)

The Bayesian data fusion approach relies on the hypothesis that one can decompose the spatial component (i.e. the prior distribution) from the other information sources (i.e. the likelihood function). Subsequently, one can assume that, conditionally to the true underlying variable, the different information sources are independent. In the Bayesian framework, this implies directly that the likelihood function can be decomposed in a product of conditional distributions. The Bayes’s theorem can thus be applied a second time in order to express the posterior distribution of a given variable Z₀ given the other secondary variables Y_i either as a function of the conditional distribution or. It is beyond the scope of this paper to present the whole underlying theory of BDF. For more details about the general description of the approach, the reader can consult [15,18] for a specific implementation of the approach.

The general BDF equations can be simplified to simple analytical formulas when it is assumed that the distributions of errors obtained by kriging and by the statistical model are Gaussian [18]. The final predicted mean nitrate concentration μ_P and its variance are then given by

(3)

where μ_k, μ_m and μ₀ are the means associated with the kriging prediction, the statistical model prediction and a “rough” estimation of the local mean obtained from the inverse distance method, respectively, and where, and are the variances associated with the kriging prediction (defined as the variance of prediction Var(Z_est,0 – Z₀), where Z_est,0 is the predicted value and Z₀ the observed value), the variance associated with the statistical model prediction (defined as the variance of the data at the end of each branch of the regression tree) and the sill of the semi-variogram, respectively.

3.4. Comparison of Methods

The presented interpolation methods yield different cartographic results because of their inherent characteristics. Each map is subject to prediction errors. To elucidate those errors and to illustrate the estimation accuracy, a “leave-one-out” cross-validation approach was performed. The following indicators were chosen for accuracy assessment:

The root mean squared error (RMSE):

(4)

The mean absolute error (MAE):

(5)

The mean error (ME):

(6)

The coefficient of determination (R²):

(7)

where n is the number of observations, is the mean value across all observations y_i and ŷ_i are the predicted values.

It is worth noting that these quality indicators only reflect the general regional accuracy of the predictions. Local improvements of the mapping process must therefore be observed directly on the maps of uncertainties.

All calculations were performed in Matlab^TM R2008b, using the geostatistical BMElib package [5]. Data preparation was carried out in ArcGIS^TM 9.3.

4. Results

4.1. Kriging

The first presented mapping method relies on the interpolation of the measured nitrate concentrations through ordinary kriging. An experimental semi-variogram was modeled as the sum of a nugget effect of 67.2 mg²/L² and an exponential model with total theoretical variance and range of 332.9 mg²/L² and 12090.8 m, respectively (Figure 1). About 66.5% of the study area has kriged nitrate concentrations lower than the standard of 50 mg/L, while nitrate concentrations exceeding 80 mg/L are observed in some parts of the study area (Figure 2(a)). The predicttion errors of the kriged concentration is minimum at the location of the monitoring stations and raises up to 20 mg/L were the density of the monitoring stations is the lowest (Figure 2(b)).

Kriged nitrate concentrations were compared to measured nitrate concentrations by a “leave-one-out” cross validation (Figure 3). The data roughly scatter around the 1:1 line. This bad prediction quality is partly due to

the very low density of the sampling network in some parts of the study area. As a consequence, when the “leave-one-out” cross validation process removes these isolated points from the dataset, the prediction solely relies on data located far away from that point, thus resulting in an inaccurate prediction and a global poor coefficient of determination. Furthermore, the presence of a nugget effect of about 20% of the total variance is also responsible for poor prediction results.

4.2. Regression Tree

4.2.1. Development of the Regression Tree

Nitrate concentrations measured at the 99 monitoring stations in January and February 2009 were related to the 13 variables listed in Table 2 through a regression tree model (Figure 4(b)). As shown by [28] such models show highly complex interaction patterns, suggesting that it is a complex combination of variables that explains ob-served pollution levels, rather than single explicative variables.

Figure 4(a), which represents a simplified version of the regression tree with only the upper part of the tree,