Jakob J Verbeek

Nikos Vlassis

Gaussian fields for semi-supervised regression and correspondence learning

Keywords: Gaussian fields, regression, active learning, model selection

Gaussian fields (GF) have recently received considerable attention for dimension reduction and semi-supervised classification. In this paper we show how the GF framework can be used for semi-supervised regression on high-dimensional data. We propose an active learning strategy based on entropy minimization and a maximum likelihood model selection method. Furthermore, we show how a recent generalization of the LLE algorithm for correspondence learning can be cast into the GF framework, which obviates the need to choose a representation dimensionality.

Introduction

Learning from high dimensional data is a challenging problem that becomes increasingly important as large and high dimensional data collections need to be analyzed in different application domains. Such data collections appear in home and industrial settings, for example in document collections, image databases, and audio recordings. In this paper we consider Gaussian fields (GFs) defined over high dimensional input data for two learning problems.

The first problem is semi-supervised regression on high dimensional data that exhibits a low dimensional structure. For example, the input data can be images of a face looking in different directions and we want to map these images to pose parameters. If we consider the images as a set of pixel values forming a vector, then the input vectors have a dimensionality that equals the number of pixels in the image. However, in this high dimensional space the data are confined to an embedding of a low dimensional manifold that contains the images of the face looking in different directions. In this setting, unsupervised data (with unknown pose parameters) can be exploited to identify the low dimensional manifold structure of the data. Due to the low dimensionality of the manifold, we expect that only few supervised examples are needed to successfully learn a function on the manifold to predict the pose parameters.

The second problem we consider is correspondence learning. In this problem two high-dimensional data sets are given that share the same underlying modes of variability. For example, the sets can be collections of images of different views of objects x and y, or the sets can be collections of images of persons x and y showing different facial expressions. The manifolds on which the two data sets lie, parameterized by the shared modes of variability (camera pose or facial expression), can be aligned if some 'correspondences' are given. A correspondence is a pair that consists of a point in the first data set and one in the second, which are known to share the same underlying low-dimensional coordinate (e.g. images of the first and second object viewed from the same direction, see Fig. 4). The goal is to predict for a data point in one set the corresponding observation in the other set. This problem is similar to regression, except that in this case the 'output' is also high dimensional. Also in this setting, the low dimensional manifolds can be recovered using unsupervised data, and can then be aligned using only relatively few correspondences due to the low dimensionality of the manifolds.

In order to successfully learn regression functions from high dimensional data either an estimation bias must be introduced by using a limited class of regression models (e.g. linear models), or alternatively, large quantities of data are required to reduce the estimation variance induced by larger model class. Recently, non-parametric techniques have been proposed [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF][START_REF] Joachims | Transductive learning via spectral graph partitioning[END_REF][START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF][START_REF] Belkin | Semi-supervised learning on Riemannian manifolds[END_REF] for unsupervised dimension reduction and semi-supervised classification that are based on a nearest neighbor graph on the set of high-dimensional 'inputs'. The graph on the data encodes the assumption that nearest neighbors in the input space tend to have similar 'outputs' (low-dimensional coordinates for the inputs in case of dimension reduction, and class labels in the case of semi-supervised classification). These techniques are based on a non-negative energy function that sums terms that are quadratic in the outputs of neighbors in the graph. With the energy function we can associate a Gaussian field (GF); a probability density function on the outputs proportional to the exponent of the negative energy function. The associated GF makes it possible to quantify uncertainty in the predictions which is useful for active learning. The main difference between the graph-based approaches and more traditional approaches lies in how the learning bias is introduced. Graph based methods are biased towards smooth functions over the graph; in this manner the bias is derived from -and reflects the structure of-the given inputs. Traditional techniques introduce a bias by limiting attention to some parametric model class. Notably, the choice of a specific function class is often motivated by reasons of computational efficiency rather than the suitability of this class for the given problem.

The rest of this paper is organized as follows. In the next section we review GFs defined on nearest neighbor graphs. In Section 3 we discuss related methods for unsupervised dimension reduction and semi-supervised classification. In section Section 4 we focus on the first learning problem: semisupervised regression using GFs. We propose a maximum likelihood model selection method to set the parameters of the GF. Furthermore, we propose an entropy minimization based query selection procedure for active learning. We present experimental results that show the effectiveness of both procedures and the advantage of the GF method over a method that separates unsupervised dimension reduction and supervised learning. In Section 5 we turn to correspondence learning using GFs. We propose a new approach for the correspondence learning based on GFs, and experimentally compare it to an existing approach [START_REF] Ham | Semisupervised alignment of manifolds[END_REF] based on the Locally Linear Embedding (LLE) algorithm [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF]. The GF approach performs as good as or better than the LLE-based method, while the GF-based approach has one parameter less: it does not require the setting of a representation dimensionality. We end the paper with conclusions and a discussion in Section 6.

Gaussian fields defined over nearest neighbor graphs

In this section we review Gaussian fields defined over neighborhood graphs. Given a set of high-dimensional 'inputs' X = {x i } n i=1 , the first step in defining the GF is to find the k nearest neighbors of each x i . Throughout this paper we define nearest neighbors on the basis of Euclidean distance (selecting the value of k is the topic of Section 4.1). Note that it is not necessary to exhaustively compute all pairwise distances to find the nearest neighbors, see e.g. [START_REF] Karger | Finding nearest neighbors in growth-restricted metrics[END_REF].

Let y i denote some scalar 'output' of x i , and let y = (y 1 , . . . , y n ) ⊤ denote the n-dimensional vector containing all outputs. Given the nearest neighbors, we define an energy function E(y) by summing over all input pairs x i and x j which are nearest neighbors, the squared difference of their outputs (y iy j ) 2 times a weighting factor a ij . If we let a ij = 0 only for nearest neighbors, then:

E(y) = n i=1 n j=1 a ij (y i -y j ) 2 = y ⊤ (D -A)y = y ⊤ Ly, (1)
where the degree matrix D is a diagonal matrix containing the row sums of the adjacency matrix A. The matrix L is known as the graph Laplacian of A. Note that if ∀ i,j a ij ≥ 0 then ∀ y E(y) ≥ 0, i.e. L is positive definite. However the converse is not true: matrices A with negative entries exist with ∀ y E(y) ≥ 0. For example, in the Locally Linear Embedding (LLE) algorithm [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] the following adjacency matrix is used

A = W + W ⊤ -W ⊤ W, (2)
where each row of W sums to 1 and w ij = 0 only for nearest neighbors. Using this adjacency matrix we have D = I, and the energy E(y) can be written as:

E(y) = n i=1 (y i - n j=1 w ij y j ) 2 ≥ 0. (3)
Clearly, the corresponding Laplacian L is positive definite, while some entries of A can be negative (typically for pairs that are not neighbors in the graph but share one or more common neighbors).

The GF that we will define over the outputs is obtained by taking the exponent of the negative energy function. In this way the a zero mean Gaussian density is obtained, whose inverse covariance matrix is given by the coefficient matrix of the energy function. However, note that the energy functions defined in ( 1) and ( 3) do not define a proper GF, since L is singular (the constant vector 1 is an eigenvector with eigenvalue zero). By adding a regularizer we ensure positive definiteness:

E(y) = y ⊤ (L + αI)y = y ⊤ My. (4)
It is easy to see that M has the same eigenvectors as L, with the eigenvalues increased by α. In effect, the regularizer replaces the infinite variance in the direction 1 (if L is used as inverse covariance matrix) of the GF density with variance α -1 . Thus, for the given points X we define the GF density over outputs y based on the energy (4) as:

p(y) = N (y; 0, β -1 C) ∝ exp - β 2 E(y), (5)
where β is a scale parameter that controls the smoothness of the density and C denotes the inverse of M. Note that although the inverse covariance matrix M is sparse, it generally induces a full covariance matrix β -1 C.

The joint GF density over all outputs induces conditional and marginal distributions over sets of outputs, which will be used for active semi-supervised learning in Section 4. In semi-supervised learning the value of outputs y i of some data points x i are known and thus fixed to these values. Without loss of generality we can assume that all indices of supervised points are smaller than those of the unsupervised points. We can thus use the partition y = [y ⊤ s y ⊤ u ] ⊤ , where y u denotes the sub-vector for unsupervised outputs and y s the subvector for the supervised outputs. Let M ss , M us and M uu denote the corresponding blocks of M:

M =   M ss M su M us M uu   . (6)
Throughout we will use this subscripting to refer to blocks of matrices and vectors that correspond to supervised and unsupervised points. Using the partitions, the energy (4) can be expanded as:

E(y) = y ⊤ u M uu y u + y ⊤ s M ss y s + 2y ⊤ u M us y s , (7)
from which it is easy to find the conditional density

p(y u |y s ) ∝ exp - β 2 y ⊤ u M uu y u + 2y ⊤ u M us y s . (8)
This conditional density can be used to predict the unsupervised outputs by computing the mean y * u (which equals the mode) of this density:

y * u = -M -1 uu M us y s . (9)
We note here that to solve [START_REF] Mackay | Bayesian Methods for Adaptive Models[END_REF] in practice, we can use the (sparse) Cholesky decomposition M uu = R ⊤ R, rather than inverting M uu (which yields a full matrix). We then first solve for R ⊤ z = -M us y s , and then for Ry * u = z. Several quantities which need to be evaluated in the techniques developed in the next sections can be computed using Cholesky factors in a similar manner.

GFs for dimension reduction and classification

In this section we discuss related non-parametric methods for dimension reduction and semi-supervised classification.

Two recent approaches for unsupervised non-linear dimension reduction can be cast in the GF framework. The LLE algorithm [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] and the Laplacian Eigenmap algorithm [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF] perform dimension reduction by minimizing E(y), as defined in (3) and (1) respectively. Note that there is a trivial solution y = c1, which is an eigenvector of L with eigenvalue zero. In order to prevent the trivial solution we can minimize the energy under the constraint that the entries in y are zero mean and have unit variance. The minimizer is then given by the eigenvector of L with the second smallest eigenvalue. To extract (d + 1) rather than d features, one can minimize E(y) under the constraint that y d+1 should be orthogonal to the already extracted eigenvectors {y i } d i=1 . The resulting minimizing vectors y i are all eigenvectors of L, with increasing eigenvalues.

Given the relationship between L and the covariance matrix of the GF, these algorithms can thus be interpreted as computing the principal components of the GF density over the outputs. Within the space of all possible output vectors y, the principal components span the linear subspace which contains most of the output variation under the GF density.

For two-class semi-supervised classification, the value of the known outputs in y s is either 0 or 1, denoting the class membership. In [START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF], classification is performed by minimizing the energy E(y), as defined in [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF], with respect to y u for fixed y s . The minimizer y * u is characterized by the linear equation:

L uu y * u = -L us y s . (10)
The unlabelled points are then assigned to class 0 if their value in y * u is smaller than 1/2, otherwise they are assigned to class 1. If we ignore the regularization with α, this approach computes the mean of the conditional density p(y u |y s ) and thresholds the values of the conditional mean at 1/2.

In other work on semi-supervised classification [START_REF] Belkin | Semi-supervised learning on Riemannian manifolds[END_REF], y is constrained to be in the span of the m eigenvectors of L with smallest eigenvalues, i.e. y = Vw where V contains the m eigenvectors as columns. The squared error y s -V s w 2 is minimized w.r.t. w, where V s contains the rows of V corresponding to the supervised points; let w * denote the minimizer. Classification is performed by computing y * u = V u w * , the values in y u are again thresholded at 1/2. In effect this approach first performs unsupervised non-linear dimension reduction using all data, as in [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF]. Second, in the low dimensional representation obtained in the first step, a linear classifier is found from the supervised data.

In [START_REF] Joachims | Transductive learning via spectral graph partitioning[END_REF] the supervised labels are not used as a hard constraint. Rather, if we let γ s denote the supervised labels, the following error function is minimized:

1 E J (y) = y ⊤ Ly + (y s -γ s ) ⊤ Γ(y s -γ s ). (11)
The first term encourages solutions that are smooth in the sense of having similar values for nearest neighbors in the graph, while the second term favors solutions that comply with the supervised labels. The trade-off between smoothness and agreement with the supervised labels, as well as different misclassification costs for each supervised data point, can be encoded in the diagonal matrix Γ.

Active semi-supervised regression with Gaussian fields

In this section we present two contributions in using GFs for semi-supervised regression. First, we propose a MAP criterion to automatically set the model parameters β and k. Second, we propose a minimum entropy query selection procedure for active learning. For clarity, we focus on regression using GFs for a single response variable. The presented material is straightforwardly generalized to prediction multiple response variables independently using the same Gaussian field.

In the related work on semi-supervised classification that we reviewed in Section 3 [START_REF] Joachims | Transductive learning via spectral graph partitioning[END_REF][START_REF] Zhu | Semi-supervised learning using Gaussian fields and harmonic functions[END_REF][START_REF] Belkin | Semi-supervised learning on Riemannian manifolds[END_REF], the adjacency matrix A that determines the energy function was specified directly by assigning non-negative entries for nearest neighbors. As we argue next, for regression it is more useful to define A through a weight matrix W, as in [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF]. We set w ij = 1/k if x j is among the k nearest neighbors of x i and w ij = 0 otherwise. In this way the energy function sums for each point the squared difference of its output and the average of its neighbors, cf.

(3). In comparison, if we directly set the entries A corresponding to nearest neighbors to 1/k then the energy function sums for each point the sum of squared differences between its output and that of its neighbors.

The apparently small difference to define A has an important consequence for extrapolation performance of regression using the GF. The effect is best understood through a simple example illustrated in Fig. 1. The input data consists of 300 data points along a spiral. The desired output of each point is its position if the data points are numbered along the spiral. For three data points (shown as circles) the desired output is specified. When A is specified directly, there is no force that drives extrapolation: if all points toward the end of the spiral are assigned the same output as the last supervised point, the contribution to the energy induced by these points is zero. In contrast, when A is defined in the LLE-based manner c.f. ( 2), the data points towards the end of the spiral are assigned higher values to ensure that the average at the last supervised point best matches the supervised value. This is necessary since the unsupervised points at the other side of the last supervised point are assigned smaller values for the interpolation between the supervised points.

Because of the lack of extrapolation in the direct method, we use the LLE based energy function [START_REF] Belkin | Laplacian eigenmaps for dimensionality reduction and data representation[END_REF] in the remainder of this paper, with w ij = 1/k when x j is one of the k neighbors of x i . Note that the extrapolation effect does not appear in the case of semi-supervised classification, since then desired outputs are always in the range of the supervised outputs (i.e. between zero and one).

Model selection

The GF density ( 5) is parameterized by k, the number of nearest neighbors connected to each point, and the variance scale parameter β. Since it is hard to set these parameters by hand, it is desirable to find optimal values automatically. We can use the maximum a-posteriori (MAP) criterion to select appropriate values for k and β.

The posterior distribution on the parameters given the supervised outputs collected in y s is p(k, β|y s ) ∝ p(y s |k, β)p(k, β). Below we assume a uniform prior p(k, β) over the parameters. In this case MAP parameter estimation reduces to maximum marginal likelihood parameter estimation, i.e. find arg max k,β p(y s |k, β). Note that if the weights w ij are not fixed to 1/k but set to maximize the marginal likelihood, then a uniform prior on the parameters cannot be used. This is because the set of models with k + 1 neighbors then includes the set of models with k neighbors as a subset.

Since p(y|β, k) is Gaussian, the marginal likelihood of the supervised points is obtained relatively easily. Using n s to denote the number of supervised points, the marginal log-likelihood of y s , given k and β, is given by:

log p(y s |β, k) = - 1 2 log |C ss | -n s log β + βy ⊤ s C -1 ss y s , , (12)
where C ss depends on k. Through differentiation, the maximum of the marginal log-likelihood over β and the maximizer β * can be analytically obtained as:

β * = arg max β log p(y s |β, k) = n s /(y ⊤ s C -1 ss y s ), (13)
which can be substituted in [START_REF] Verbeek | Non-linear CCA and PCA by alignment of local models[END_REF] to find

l * = log p(y s |β * , k) = - 1 2 log |C ss | + n s + n s log(y ⊤ s C -1 ss y s /n s ) . (14)
Thus, to find the optimal number of neighbors k * we can evaluate (14) for k = {1, 2, . . . , k max } and pick the value of k for which ( 14) is maximum.

To evaluate (14), we need to find C ss , which can be computed without explicitly computing all entries of C. Recall that C = M -1 , then using the Cholesky decomposition M = R ⊤ R we can find the required entries of C ss by noting that

C ij = e ⊤ i Ce j = e ⊤ i (R ⊤ R) -1 e j s ⊤ i s j
, where e i is a vector with all zeros except for the i-th entry which equals 1, and s i = R -⊤ e i . Since R is triangular we can solve for s i in time linear in n.

Maximization of the marginal likelihood over the regularizer α cannot be performed in closed form, but it is possible to use gradient based techniques. The derivative of l * (the marginal log-likelihood maximized over β) with respect to the regularizer α equals:

∂l * ∂α = 1 2 Tr{C -1 ss [C 2 ] ss } - n s 2 y ⊤ s C -1 ss [C 2 ] ss C -1 ss y s y ⊤ s C -1 ss y s . (15)
To find the required block of C 2 it is not necessary to explicitly compute this product. To see this, note that [C 2 ] ij = e ⊤ i CCe j s ⊤ i s j . Again using the Cholesky decomposition of M (ignoring β for clarity), we can write s i = Ce i = (R ⊤ R) -1 e i . We can solve for s i by first solving R ⊤ z = e i and then for Rs i = z. However, after each gradient step the Cholesky factors R will need to be recomputed. This makes optimization over the regularizer computationally demanding. In our experiments we therefore fixed the regularizer at α = 10 -11 .

Query selection

The amount of supervised data available in practice is often limited, since it is generally costly to acquire. This motivates the use of active learning approaches, which determine, on the basis of a set of unsupervised and supervised data, for which unsupervised data the supervised signal is most informative. In this way the user's effort is not wasted on answering uninformative queries, as could be the case if supervision is queried for a randomly selected inputs.

To determine which query set of fixed size is the most informative, we may consider the entropy H(y u |y s ) of the conditional density of the unlabelled points given a labelled query set, and aim for the query set that leads to the lowest conditional entropy (see also [START_REF] Mackay | Bayesian Methods for Adaptive Models[END_REF][START_REF] Hwa | Sample selection for statistical grammar induction[END_REF]). Intuitively, we want the uncertainty of our predictions on the unlabelled points, given the labelled points, to be as small as possible. The rationale is that, if we assume that the mean of the field will approach the true value of the response variable, then the distribution p(y u |y s ) should be concentrated as tight as possible around its mean. The chain-rule of entropies tells us that H(y) = H(y s ) + H(y u |y s ). Since the entropy H(y) of the complete field is fixed, we can equivalently aim for the query set with the largest entropy. To compute the entropy H(y s ) = 1 2 log |C ss |-1 2 n s log β+const. we use the Cholesky factors of M as above to find C ss .

To find the single variable with maximum entropy (i.e. maximum variance) we should compute all marginal variances: the diagonal elements of C. However, already by selecting the maximum variance variable from a random subset of all variables, we have large probability of selecting one of the variables with largest variance. In fact, as in [START_REF] Smola | Sparse greedy matrix approximation for machine learning[END_REF], we can compute confidence bounds of the type: With probability at least 1δ the variable with largest variance among m uniformly randomly selected variables is among the ε fraction of variables with largest variance, if m ≥ ⌈log δ/ log (1ε)⌉. For example, if we use a random subset of 59 variables then with probability at least 95% the selected variable will be among the 5% of variables with the largest variances.

The number of possible query sets grows quickly with n s and heuristic search methods are inevitable already when n s is moderately large. Given some initial query set, a simple, and in practice very successful, scheme to improve the query set is the following. Until some stopping criterion is met, randomly some variable y i not in the current query set is selected. An exchange is made if replacing one of the variables in the current query set with y i yields a higher joint entropy of the query set. To generate the initial query set either a random subset can be used or a greedy procedure, that adds new elements to the query set one-by-one, which we describe below. When the size of the query set is not known in advance, it is also possible to use the greedy procedure by itself. In that case, when the user decides to label another data point, the query that optimally augments the current query set is determined.

Given an initial query set s ⊂ {1, . . . , N } we want to find the best query i to add to s, with i ∈ u = {1, . . . , N }\s such that the entropy in the remaining variables H(y u\i |y {i}∪s ) = H(y u |y s ) -H(y i |y s ) is minimized. To do this we find the i ∈ u for which H(y i |y s ) is maximum; to compute these entropies we use the Cholesky decomposition of M uu . Note that this decomposition is also used to predict y u from y s . In cases where the Cholesky decomposition of M uu is not available (because y u is not predicted or predicted without use of the Cholesky decomposition), the decomposition of M can used instead. To see this, note that from the chain-rule of entropies we have that H(y u\i |y {i}∪s ) = H(y)-H(y {i}∪s ). Thus to find the optimal query to add to s, we find the query i that, together with s, has the maximum entropy H(y {i}∪s ) in the complete field. The entropies H(y {i}∪s ) are computed from the Cholesky factors of M.

Experimental results

First we consider a face-pose estimation experiment to give a qualitative illustration of the techniques described above. The input data consists of a set of 2000 images (40×40 pixels) of a face that looks in different directions, varying both horizontally and vertically (see the images depicted in Fig. 2). 2 The goal is to predict the pose of unsupervised images given only 10 supervised images.

A neighborhood graph with 20 neighbors was used, based on Euclidean distance. The supervised data was obtained as follows: first the user has to place a random image in a two-dimensional plane. Then, using the minimum entropy principle a next query image is selected and the user places this image also in the plane. This process was repeated until ten images were placed in the plane. The two coordinates of the position in the plane of each image are considered as the supervised values of the two response variables which correspond to the horizontal and vertical viewing direction (pose) of the face.

In Fig. 2 the predicted pose parameters of all images are shown, while the coordinates of the 10 supervised images are indicated with the digits 1, . . . , 10. Observe that the query selection procedure mainly requested supervision for extreme poses. To illustrate the predictive power of the model, we sampled ten equally spaced points (depicted as circles) along four lines in the pose space. The images nearest to the sampled poses are plotted next to the figure. These images show that using only 10 supervised examples it is possible to accurately predict the pose depicted in the remaining unsupervised images.

We continue with quantitative experimental results obtained using a data set of 698 images (64×64 pixels) of a face rendered by computer graphics engine. 3The face is viewed and lit from different directions and for all images the ground truth viewing direction (longitude and latitude) and lighting direction (longitude) are available. The task is to predict the viewing direction and lighting direction of unsupervised images. The three response variables were normalized to the range [0, 3.5] yielding approximately unit variance.

In the experiments we used all 698 images, some of which were supervised. The supervised images were selected either randomly or by an active learning scheme. The squared errors between the predictions and the ground truth are averaged over the unsupervised points and over the different response variables. All reported results are averages over 20 experiments and we used t-tests to assess the p-value significance levels of the obtained results presented in Fig. 3 (note the log scale).

Query selection

In Fig. 3 we show both for random query selection (solid lines) and for active query selection (dotted lines) the squared errors against the number of neighbors k used in the model. The different lines (from top to bottom) show results for 10, 20 and 100 labelled points. We started the query selection procedure with the greedy query selection method, i.e. selecting one-by-one the queries that are expected to yield the most information given the other already selected query points. After the greedy initialization we improved the query set by an iterative procedure. We randomly selected one point not in the query set and computed whether the query set entropy increased if this point replaced any of the points in the current query set. We terminated the iterations if no entropy increase was achieved after considering 20 unsupervised points. 3. Errors on unsupervised points against number of neighbors, using 10, 20 and 100 randomly (solid lines) and actively (dotted lines) selected queries (top panel). Some example images (bottom right). Model selection results versus using optimal k value (table ).

We ran experiments for neighborhood graphs of different connectivity: k = 2, 3, . . . , 20. In all but 4 of the 57 combinations of k and a number of labelled points the active learning algorithm performs significantly (p < 0.001) better than using random queries. Also note that the standard deviation in the performance of the active learning approach is much smaller. We conclude that the active learning algorithm outperformed the usage of random queries.

Model selection

Next, we consider the performance of the maximum likelihood model selection procedure described in Section 4.1. Note that query selection depends on the number of neighbors that are used, and that model selection depends on the supervised queries. To combine model selection with query selection we use a greedy procedure: first a random query is labelled, and then alternatingly k is selected using the labelled points, and a new query is selected depending on the current k and already labelled points. If random queries are selected, we can perform model selection directly using all labelled points.

The table in Fig. 3 shows the average results that were obtained with the optimal value of k in each of the 20 experiments and with the value of k found by model selection. In each of the 20 experiments the optimal value of k was determined separately when combined with active and random queries: it is the number of neighbors for which the smallest prediction error was obtained (using active or random queries respectively).

Using random queries, model selection is near optimal: the optimal k does not give significantly better results than using the automatically selected k (p = 0.1). However, for actively selected queries and only 10 supervised points, the results for optimal k are significantly better (p < 0.02). The reason is that for few supervised points, model selection picks too small k which causes the query selection procedure to perform poorly. Comparing active and random queries, we found that for 10 queries there is no significant difference between active and random selection if model selection is used (p = 0.4). For more queries active selection is significantly better (p < 10 -4 ).

Comparison to staged unsupervised and supervised learning

In this experiment we compare the GF approach to the staged learning approach of [START_REF] Belkin | Semi-supervised learning on Riemannian manifolds[END_REF] for semi-supervised classification as discussed in Section 3. The same data set as in the previous experiment was used, and each experiment was repeated 20 times with supervised points chosen uniformly random 4 in each experiment. We computed the prediction error on the remaining unsupervised points for different numbers of supervised points n s , and numbers of nearest neighbors k to construct the graph. For the staged learning approach we have to set a number of eigenvectors m that will be used. In Tab. 1 we report the results for the optimal number of eigenvectors m ∈ {1, . . . , 20}, i.e. the number that led to minimum average error. The results are printed bold if there is a significant difference in performance. The performance difference between Prediction errors ×10 (average and st.dev.) obtained using GF and staged learning (SL). The optimal number of eigenvectors for SL is also given for each setting.

the two methods is mostly not significant (p < 0.01) for n s ∈ {10, 20}. For n s = 100 supervised points the GF method performs significantly (p < 0.01) better. Note that the GF method does not require the choice for a particular number of eigenvectors, and thus has one parameter less.

Correspondence learning

The second learning problem we consider is learning correspondences between two sets of data points. Each data set is assumed to be sampled from a different high-dimensional embedding of the same underlying low-dimensional manifold. The learning examples consist of (i) several supervised pairs of points (one from each data set) that are in correspondence, i.e. have the same (unknown) coordinate on the underlying manifold and (ii) many unsupervised points in each set for which the corresponding point in the other set is unknown. The goal is to predict other correspondences, i.e. for unsupervised points we want to predict which point in the other set corresponds to it.

Correspondence learning with the LLE algorithm

The approach to correspondence learning taken in [START_REF] Ham | Semisupervised alignment of manifolds[END_REF] is to perform dimension reduction for both sets simultaneously with the LLE algorithm, explicitly aiming at discovering the underlying low-dimensional coordinates. The energy function adds the LLE energy functions (3) of each set. As we explain next, the correspondences are accounted for by constraining the low-dimensional coordinates of corresponding points to coincide.

Without loss of generality we order the points in both sets in such a manner that the corresponding points have equal indices in each set, and all points with a correspondence have a smaller index than points without a correspondence. Let us use a and b to index the two data sets. Let L a and L b be the Laplacians, c.f. ( 1)-( 3), of the separate sets as used by LLE. We use the partitions

L a =   L a ss L a su L a us L a uu   , y a =   y a s y a u   , (16)
and similar partitions for L b and y b . The two Laplacians are combined to form

L =     L a ss + L b ss L a su L b su L a us L a uu 0 L b us 0 L b uu     . (17)
Taking into account the constraints y a s = y b s y s , and letting y = [y ⊤ s y a ⊤ u y b ⊤ u ] ⊤ , the sum of the energy functions, c.f. c.f. ( 1)-( 3), of the two sets equals:

E(y) = y ⊤ Ly. (18)
Analogous to the normal LLE algorithm, the eigenvectors of L with smallest eigenvalues give the minimum energy assignments of y a u , y b u , and y s . The low-dimensional coordinates of each point are in the corresponding row of the matrix of which the columns are the smallest eigenvectors of L. To predict a correspondence of an unsupervised point in set a, the nearest (using Euclidean distance in the low-dimensional representation) point from set b is found. Note that the dimensionality of the low-dimensional space, as well as relative scaling of the dimensions, influences the obtained results. In some cases an appropriate dimensionality for the embedding is known. In other cases, however, this dimensionality has to be estimated. Furthermore, in the above approach an appropriate number of nearest neighbors has to be determined for each set.

The semi-parametric correspondence learning setup in [START_REF] Verbeek | Non-linear CCA and PCA by alignment of local models[END_REF] differs in that it delivers a mixture of factor analyzers from which the high-dimensional corresponding observations can be reconstructed without storing the training data. However, in the semi-parametric approach some trade-off has to be made between satisfying the correspondences (ensuring that correspondences have similar low dimensional coordinates) and mutual alignment of the factor analyzers fitted to each set of observations (ensuring that each factor analyzer assigns approximately the same low-dimensional coordinates to data points).

Correspondence learning with Gaussian fields

By casting the correspondence learning in the GF framework, we do not need an explicit representation in a low-dimensional space. Essentially, we compute quantities that are averaged over possible representations based on their likelihood under the GF distribution. As a first step we define a GF for each point set separately as in the previous sections, with inverse covariance matrices βM a and βM b . We may view this as a single GF with two independent subfields. To account for correspondences, we constrain variables that are linked by the correspondences to have identical (but unknown) values. The constraints induce correlations between other variables in the subfields. The inverse covariance matrix βM is obtained in the same way as L above, where the two variables of each correspondence pair have been replaced by a single variable. Again we use C to denote the inverse of M of the combined field. The expected squared difference (y iy j ) 2 under the constrained GF distribution is:

E[(y i -y j ) 2 ] = E[y 2 i ] + E[y 2 j ] -2E[y i y j ] = β -1 (C ii + C jj -2C ij ). ( 19
)
Given a query point i from one set we can then find the point j in the other set with smallest expected squared difference (note that this point j is invariant for the value of β).

The LLE based method could be viewed as an approximation to our GF based method in the following sense. Let v t and λ t denote eigenvectors and the corresponding eigenvalues of C, then we can write

C ij = t v t (i)v t (j)λ t . (20)
Let the eigenvalues be ordered from large to small: λ t ≥ λ t+1 , the LLE based method [START_REF] Ham | Semisupervised alignment of manifolds[END_REF] is obtained (for embedding dimensionality d) as follows. We set λ t ← 1 for t ∈ {1, . . . , d + 1} and λ t ← 0 for t > d + 1. The squared distance in the eigenvector embedding used by the LLE based method equals C ii + C jj -2C ij , where the C ij are computed as in (20) with the replaced λ t .

Model selection

In the correspondence problem we can use a marginal likelihood model selection procedure (selecting an appropriate number of neighbors k and variance scale parameter β) similar to the procedure used in the previous section. However, here we consider the high-dimensional input vectors as the variables that we want to predict. For the n s given correspondences both high-dimensional input vectors belonging to sets a and b are known.

The model selection method of the previous section is based on maximizing the marginal log-likelihood of the supervised data. Suppose that the input vectors from set a and b have dimensionality respectively d a and d b . Let d = d a + d b , then we define the n s × d matrix Z as the matrix having in the i-th row the concatenation of the input vectors of the i-th correspondence pair. In the current setting, model selection will be based on maximizing the marginal log-likelihood of Z. Assuming that the columns z j of Z are independently and identically distributed5 according to the GF distribution, this translates into maximizing the sum of the marginal log-likelihoods of the columns of Z. Thus model selection is based on finding k and β that maximize:

p(Z|k, β) = d j=1 log p(z j |k, β) (21)
= d j=1 - 1 2 log |C ss | -n s log β + βz ⊤ j C -1 ss z j (22) = - 1 2 d log |C ss | -dn s log β + βTr{Z ⊤ C -1 ss Z} . (23)
The maximization with respect to β can again be done in closed form yielding an expressions similar to ( 13) and ( 14):

β * = arg max β log p(Z|β, k) = dn s /Tr{Z ⊤ C -1 ss Z}, (24)
l * = log p(Z|β * , k) (25) = - 1 2 d log |C ss | + dn s + dn s log Tr{Z ⊤ C -1 ss Z}/(dn s ) . (26)
To find the value of k for which l * is maximal (and the corresponding β * ) we evaluate (26) for different values of k. To do this we need to find C ss , which can be computed with methods discussed in the previous section. Note that to compute the trace in (26) we can use Tr{Z ⊤ C -1 ss Z} = d j=1 z ⊤ j C -1 ss z j . Thus the required amount of computation depends only linearly on d.

Experimental results

We compare the performance of our GF approach with the LLE approach proposed in [START_REF] Ham | Semisupervised alignment of manifolds[END_REF] using representations of different dimensionality. The data sets we used each contain 2500 views of a toy puppet recorded by positioning a camera on different positions on the hemisphere centered at the object. 6 For each view of one object there is a corresponding view of the other object recorded from the same camera position, see the images in Fig. 4. Given several corresponding views, the task is to predict for each remaining image in set a which is the corresponding image in set b and vice-versa. We compare the position of the camera of the true correspondence and the predicted correspondence.

The error we measure is the summed squared error in the longitude of camera (0 • -360 • ) and latitude (0 • -90 • ). The neighborhood graph (using k = 20 neighbors) was based on Euclidean distances when between the 64 × 64 pixel images (considering each pixel as a coordinate in the image space).

In Fig. 4 the mean and standard deviation of the errors are shown of ten experiments using different random sets of 20, 50, 100 and 200 correspondences. First, we can see that for dimensionality smaller than three the LLE based method performs poorly; the differences with GF are significant with p < 0.01. The reason is that of the two degrees of freedom in the data, one (the longitude of the camera) is periodic. Therefore an additional dimension is required for a suitable representation. Second, the GF approach leads to equal or smaller errors and has less variance in the errors than the LLE based method using three or more dimensions and less than 100 correspondences (7 of the 9 differences are significant with p < 0.05). We conclude that the GF method (i) removes the need to find the optimal representation dimensionality and (ii) performs as good as or better than the LLE based method.

In Fig. 5 we show results for model selection, averaged over ten experiments, using different randomly selected sets of 50 correspondences. In each experiment we computed for k = 2, 4, 6, . . . , 20 neighbors the marginal log-likelihood l * of the 2 × 64 2 observed pixel values of the correspondences. In each experiment, we predicted the correspondences with the LLE based method and the GF method. The average prediction error and marginal log-likelihood are plotted against the number of neighbors in Fig. 5. Note that marginal likelihood gives a good indication of the prediction error, and that -on averagefor k = 8 the highest marginal log-likelihood is obtained and both methods give the smallest prediction errors. Another interesting quantity is the average prediction error when using in each experiment the value of k for which the maximum marginal log-likelihood was achieved in that experiment. This is the error that would be obtained if we used the maximum likelihood criterion for model selection. The average of this error is given by the horizontal lines in the left panel of Fig. 5. In five experiments k = 8 was selected, three times k = 10 was selected and in two cases k = 12 was selected. For both the LLE-based and the GF approach, model selection leads to prediction errors near the error obtained using the optimal value of k.

Conclusion and discussion

In this paper we have shown how sparse GFs based on a neighborhood graph can be used for semi-supervised regression and correspondence learning. For semi-supervised regression we have shown that (i) a number of neighbors that 20 yields good generalization can effectively be estimated using a maximumlikelihood criterion (ii) smaller errors can be obtained with an equal number of labelled examples if an entropy minimization principle is used to select which data points to label (iii) comparable or smaller prediction errors are obtained as compared to a staged unsupervised-supervised learning approach. For correspondence learning we have shown that our GF approach, for which we do not have to fix a representation dimensionality, performs comparable or better than the LLE based method. Furthermore, also for correspondence learning a number of neighbors that gives accurate predictions can be determined through a maximum likelihood criterion.

An alternative to using Cholesky decompositions, is to use iterative algorithms, e.g. Lanczos iterations, to solve the linear systems as required for model selection and query selection. Such methods are attractive in active learning settings since the solution after labelling a new point will be similar to the solution before labelling the point. Therefore, accurate solutions are expected to be found in relatively few iterations if the iterations are started with the (approximate) solution before the last labelling.

The Gaussian field (GF) models considered here resemble Gaussian Processes (GP) [START_REF] Neal | Bayesian Learning for Neural Networks[END_REF]. A GP is given by a mean function m(x) and a covariance function C(x, x ′ ) which together define a Gaussian density over outputs for some input domain X . The marginal Gaussian on the outputs of {x i ∈ X } n i=1 has mean vector [m(x 1 ) . . . m(x n )] ⊤ and the (i, j) th entry of the covariance matrix is given by C(x i , x j ). Thus, the data are used to specify the covariance matrix rather than the inverse covariance function as in GFs. Therefore, in a GP the marginal on the outputs of several data points does not depend on other unsupervised data points. In contrast, in a GF the marginal distribution on the outputs of several data points does depend on other unsupervised points, since the covariance matrix of the marginal is given by the corresponding block of the inverse of M and thus also depends on the parts of M corresponding to other points. In this manner unsupervised data modifies the covariance matrix to account for the manifold structure derived from the unsupervised data.

Fig. 1 .

 1 Fig.1. Input data distributed along spiral (left panel). True and predicted output for the data (right panel). The supervised points are encircled in both panels. Predictions are made using the direct and LLE based method to specify matrix A.

Fig. 2 .

 2 Fig. 2. Predicted pose parameters of 2000 images, together with the pose parameters of the ten supervised images. For 40 other points in the pose space, the image with the most similar predicted pose-parameters is shown.

  28 0.15 0.05 R S 0.31 0.16 0.05 A O 0.12 0.08 0.03 A S 0.32 0.09 0.03 R/A: random / active queries O/S: using optimal/selected k

Fig.

  Fig.3. Errors on unsupervised points against number of neighbors, using 10, 20 and 100 randomly (solid lines) and actively (dotted lines) selected queries (top panel). Some example images (bottom right). Model selection results versus using optimal k value (table).

Fig. 4 .

 4 Fig. 4. Mean (top left) and standard deviations (top right) of error against number of correspondences. Results are shown for the LLE based method using 1 up to 5 dimensions and the GF method. Three examples of corresponding images connected by bars (bottom).

Fig. 5 .

 5 Fig.5. Mean error against the number of neighbors, for LLE based method (dotted) and GF (solid). Horizontal lines give the average obtained result if model selection is used (left panel). Mean and standard deviations of marginal log-likelihood against the number of neighbors (right panel).

Table 1

 1

	ns	k	4	6	8	10	12	14	16	18	20
	10	GF	5.1	4.5	5.1	3.7	4.3	4.0	4.1	3.8	3.2
			±2.0	±2.0	±3.0	±1.3	±2.1	±2.2	±2.6	±1.7	±1.5
		SL 11.0 /1	6.3 /4	7.6 /4	3.3 /4 2.9 /4	2.8 /4	3.3 /4	4.5 /4 4.4 /4
			±1.3	±4.0	±7.5	±1.8	±1.3	±1.1	±2.4	±5.5	±10.1
	20	GF	2.7	2.2	2.1	1.9	2.0	1.8	1.8	1.8	1.8
			±0.9	±0.5	±0.6	±0.5	±0.8	±0.5	±0.6	±0.5	±0.4
		SL	3.3 /8	3.1 /7	1.8 /8	1.6 /7 1.5 /7	1.7 /7	1.8 /7	1.8 /8 1.9 /7
			±3.6	±1.2	±0.9	±0.6	±0.3	±1.1	±0.9	±0.5	±0.6
	100 GF	1.0	0.8	0.8	0.8	0.7	0.7	0.7	0.8	0.8
			±0.1	±0.1	±0.1	±0.1	±0.1	±0.1	±0.1	±0.1	±0.1
		SL 1.2 /16 1.0 /17 1.0 /16 0.9 /8 0.9 /7 0.9 /11 0.9 /11 1.0 /8 1.0 /9
			±0.2	±0.1	±0.1	±0.1	±0.0	±0.1	±0.1	±0.0	±0.1

This is a slightly simplified account; in[START_REF] Joachims | Transductive learning via spectral graph partitioning[END_REF] solutions are constrained to be linear combinations of the smallest eigenvectors of L and to satisfy y ⊤ 1 = 0 and y ⊤ y = n.

The data set is available from http://www.science.uva.nl/ ∼ jverbeek.

The data set, also used in[START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], is available from http://isomap.stanford.edu.

We use random queries to avoid biasing results due to GF-based query selection.

This assumption allows for straightforward generalization of the model selection method used in the previous section. Dropping this assumption requires non-trivial specification of conditional fields over columns of Z given other columns.

The data is available from http://ls7-www.cs.uni-dortmund.de/ ∼ peters.

Acknowledgments JJV is supported by the Technology Foundation STW (AIF4997) applied science division of NWO and the technology program of the Dutch Ministry of Economic Affairs.


A Van Geen
email: avangeen@ldeo.columbia.edu

K Radloff

Z Aziz

Z Cheng

M R Huq

K M Ahmed

B Weinman

S Goodbred

Yunlin H B Jung

Jacques Zheng

Y Zheng

M Berg

P T K Trang

L Charlet

J Metral

D Tisserand

S Guillot

S Chakraborty

A P Gajurel

B N Upreti

Comparison of arsenic

published or not. The documents may come

Introduction

Even though scientists throughout the world have been mobilized by the health impact of elevated As levels in reducing groundwater of southern Asia, attempts to explain the observed patterns from first principles have so far met with limited success. One difficulty is that it remains unclear which particulate phase(s) play a role in regulating the partitioning of As between the dissolved and particulate phase -be it through solubility control or adsorption. There is little doubt that microbial reduction of Fe oxyhydroxides is a precondition for significant release of As to groundwater [START_REF] Nickson | Arsenic poisoning of Bangladesh groundwater[END_REF][START_REF] Bgs | Arsenic contamination of groundwater in Bangladesh[END_REF]. But as pointed out by [START_REF] Polizzotto | Solid-phases and desorption processes of arsenic within Bangladesh sediments[END_REF], it is hard to reconcile the leading hypothesis linking the release of As to groundwater to the reduction of Fe(III) oxyhydroxides with the evidently very low proportion, possibly absence, of such phases even in shallow grey aquifers where dissolved As concentrations are low. Binding constants for As on Fe(III) oxhydroxides derived from controlled laboratory studies tend to underestimate the concentration of As in solution when applied to natural conditions, even if competing ligands such as phosphate and silicate are taken into account [START_REF] Bgs | Arsenic contamination of groundwater in Bangladesh[END_REF][START_REF] Swartz | Mobility of arsenic in a Bangladesh aquifer: Inferences from geochemical profiles, leaching data, and mineralogical characterization[END_REF]. In the case of a particularly well characterized shallow aquifer in Vietnam and relatively simple hydrology, [START_REF] Postma | Arsenic in groundwater of the Red River floodplain, Vietnam: Controlling geochemical processes and reactive transport modeling[END_REF] had to allow for two different stability constants for Fe oxyhydroxides, an adjustable As/Fe ratio in dissolving Fe oxhydroxides, and a reduced density of adsorption sites for Fe (III) oxhydroxides in order to model their data. This suggests that quantitatively relating laboratory-based studies of As behavior in model systems to redox changes in the natural environment, while clearly the ultimate goal, may remain difficult to achieve for a while to come.

Over the past few years, contributors to the present study have taken a more empirical approach by trying to determine if there is a systematic relationship between As concentrations in the groundwater and a plausible reservoir of exchangeable As in aquifer particles [START_REF] Swartz | Mobility of arsenic in a Bangladesh aquifer: Inferences from geochemical profiles, leaching data, and mineralogical characterization[END_REF][START_REF] Zheng | Geochemical and hydrogeological contrasts between shallow and deeper aquifers in two villages of Araihazar, Bangladesh: Implications for deeper aquifers as drinking water sources[END_REF]. A key tool developed with this goal in mind has been the needlesampler, a simple device that collects groundwater as well as sediment as a slurry from precisely the same depth interval [START_REF] Van Geen | Testing groundwater for arsenic in Bangladesh before installing a well[END_REF][START_REF] Van Geen | A transect of groundwater and sediment properties in Araihazar, Bangladesh: Further evidence of decoupling between As and Fe mobilization[END_REF]. Deploying the needle-sampler overcomes the problem when collecting groundwater from a well and aquifer sediment by coring nearby that the two intervals may be geochemically distinct because of the highly variable nature of fluvio-deltaic deposits. Using the needle-sampler in a broad range of settings in Bangladesh, [START_REF] Van Geen | Flushing history as a hydrogeological control on the regional distribution of arsenic in shallow groundwater of the Bengal Basin[END_REF] recently documented that a roughly linear relationship between dissolved As in groundwater and P-extractable As in the particulate phase held across 3 orders of magnitude in concentrations and interpreted this finding as an indication of adsorptive equilibrium. The importance of reductive dissolution as a pre-condition for As release to groundwater was acknowledged, as this relationship was shown to hold only for intervals that were sufficiently reducing.

The new evidence for adsorptive equilibrium may guide future studies of the specific mechanisms of As mobilization but, perhaps more significantly, provided the first quantitative estimate of the likely time scale over which an aquifer can be flushed of its initial As content by groundwater flow (BGS/DPHE, 2001). The distribution coefficient inferred from the data in Bangladesh (K d = 4 mL/g) shows that, on average, the As content of groundwater accounts for ~3% of the combined pools in the dissolved and mobilizable particulate phase. Converted to a retardation factor, the estimated distribution coefficient indicates that As is transported through an aquifer ~40 times more slowly than groundwater because of adsorption. Using this information and a simple advection-diffusion model, [START_REF] Van Geen | Flushing history as a hydrogeological control on the regional distribution of arsenic in shallow groundwater of the Bengal Basin[END_REF] showed that shallow aquifers deposited during the late Holocene could be flushed of their initial As content in a few hundred to a few thousand years, depending on the velocity of groundwater flow [START_REF] Bgs | Arsenic contamination of groundwater in Bangladesh[END_REF][START_REF] Mcarthur | Natural organic matter in sedimentary basins and its relation arsenic in anoxic ground water: the example of West Bengal and its worldwide implications[END_REF][START_REF] Ravenscroft | Arsenic in groundwater of the Bengal Basin, Bangladesh: Distribution, field relations, and hydrogeological setting[END_REF][START_REF] Klump | Groundwater dynamics and arsenic mobilization in Bangladesh assessed using noble gases and tritium[END_REF]. An argument was presented, based on a representative range of groundwater flow velocities, that the regional and local distribution of As in shallow aquifers of the Bengal Basin could therefore reflect different flushing histories rather than a difference in the initial As content of aquifer particles [START_REF] Meharg | Codeposition of organic carbon and arsenic in Bengal Delta aquifers[END_REF].

It is worth emphasizing that a subset of intervals sampled in Bangladesh characterized by leachable Fe(II)/Fe ratios <0.5 in the particulate phase or dissolved Fe <0.2 mg/L in groundwater containing <10 μg/L dissolved As did not fit the relationship corresponding to a K d = 4 mL/g. The outliers, characterized by a significant deficit in dissolved As compared to the prediction from the relationship and the level of P-extracted As, were tentatively attributed to the persistence of binding sites with a higher affinity for As under insufficiently reducing conditions, possibly a very small proportion of Fe(III) oxyhydroxides. Much less reducing conditions corresponding to aquifers containing orange-colored sands of likely Pleistocene age and particularly low leachable Fe(II)/Fe ratios were not considered in this analysis (BGS/ DPHE, 2001; [START_REF] Horneman | Decoupling of As and Fe release to Bangladesh groundwater under reducing conditions. Part I: Evidence from sediment profiles[END_REF][START_REF] Zheng | Geochemical and hydrogeological contrasts between shallow and deeper aquifers in two villages of Araihazar, Bangladesh: Implications for deeper aquifers as drinking water sources[END_REF].

A more recent data set obtained with the needle-sampler from an As-affected region of West-Bengal, India, has shown that the relationship between dissolved As and P-extractable As can be more complex in other settings [START_REF] Métral | Comparison of dissolved and particulate arsenic distributions in shallow aquifers of Chakdaha, India, and Araihazar, Bangladesh[END_REF]. This was attributed to generally less reducing conditions in grey shallow aquifers bordering the Hooghly River and, more tentatively, to mechanisms of As mobilization from freshly deposited sediment possibly distinct from the reductive dissolution of Fe oxyhydroxides such as the oxidation of As-rich Fe sulfides [START_REF] Polizzotto | Solid-phases and desorption processes of arsenic within Bangladesh sediments[END_REF]. In the present study, the relationship between dissolved and P-extractable As was further tested by reporting new needle-sampler data from As-affected areas in Vietnam and Nepal. On the basis of these new observations, the sensitivity of the previously proposed flushing model to a wider range of distribution coefficients is also documented.

Study sites

The geomorphological setting of 3 pairs of needle-sampler profiles collected in Bangladesh and one transect composed of 7 profiles from West Bengal, India, has been described elsewhere [START_REF] Van Geen | A transect of groundwater and sediment properties in Araihazar, Bangladesh: Further evidence of decoupling between As and Fe mobilization[END_REF][START_REF] Van Geen | Flushing history as a hydrogeological control on the regional distribution of arsenic in shallow groundwater of the Bengal Basin[END_REF][START_REF] Métral | Comparison of dissolved and particulate arsenic distributions in shallow aquifers of Chakdaha, India, and Araihazar, Bangladesh[END_REF]. In brief, an alluvial fan and an estuarine delta-plain were sampled to a depth of ~30 m in northwestern and southern Bangladesh, respectively (Fig. 1). Profiles collected in both central Bangladesh and West Bengal instead are representative of an intermediate lowland floodplain environment, with the distinction that the supply of sediment essentially ceased a few hundred years ago in the first case due to the shift of a major stream [START_REF] Weinman | Contributions of floodplain stratigraphy and evolution to the spatial patterns of groundwater arsenic in Araihazar, Bangladesh[END_REF], while alluvium continues to be supplied at the second site along the banks of the Hooghly River to this day. The degree of reduction of Fe oxydroxides leachable from the particulate phase varied with depth and location, but aquifer sands from intervals analyzed in Bangladesh and West Bengal were all grey in color and of known or likely Holocene age.

A transect of 7 new needle-sampler profiles was collected in Vietnam in April 2006 from Van Phuc village, located in the northern portion of the Red River delta and 13 km south of Hanoi (Fig. 1, [START_REF] Berg | Arsenic contamination of groundwater and drinking water in Vietnam: a human health threat[END_REF][START_REF] Berg | Hydrological and sedimentary controls leading to arsenic contamination of groundwater in the Hanoi Area, Vietnam: The impact of iron-arsenic ratios, peat, river bank deposits, and excessive groundwater abstraction[END_REF]. The Red River originates from the Yunnan Plateau in southern China, but sediment reaching the Hanoi area is derived from the entire catchment. There is a sharp contrast in As concentrations in shallow (<30 m) groundwater of the study area, with levels measured in the laboratory <10 μg/L in the northwestern portion of the village and concentrations up to 500 μg/L in the remaining portion to the south (see Fig. 6 in [START_REF] Berg | Magnitude of arsenic pollution in the Mekong and Red River Deltas-Cambodia and Vietnam[END_REF]also Eiche et al., 2008). The discussion here is limited to grey sands of likely Holocene age and data are reported for 4 needle-sampler profiles extending to 30 m depth and spaced equally along the 0.5 km southern portion of the transect (Table 1).

The location of a new 1.4 km transect composed of 5 needle-sampler profiles collected in May 2007 in the alluvial floodplain of the Terai in Nepal, near the south-central border with India, was also chosen to maximize the range of groundwater As concentrations (Fig. 1, Table 1). A survey of existing wells based on laboratory measurements indicated that groundwater As concentrations range from <10 μg/L SW of Ramgram municipality, Nawalparasi district, to 500 μg/L to the SE (see p. 94 in [START_REF] Shrestha | The State of Arsenic in Nepal-2003, National Arsenic Steering Committee[END_REF]. Whereas aquifer sands recovered from depths up to 25 m in this area were also grey in color and probably deposited during the Holocene, their origin can be linked to the erosion of the poorly consolidated sedimentary rocks of the Siwalik Group that constitute the Churia hills to the north [START_REF] Willis | Fluvial systems in the Siwalik Miocene and Wyoming Paleogene[END_REF][START_REF] Upreti | An over view of the stratigraphy and tectonics of the Nepal Himalaya. Geology of the Nepal Himalayas: Recent Advances[END_REF][START_REF] Guillot | Bengal arsenic, an archive of Himalaya orogeny and paleohydrology[END_REF]. The transect is too distant from any of the major rivers draining the Terai to receive any direct of sediment input today from the Lesser, Higher or Tethys Himalaya that lie to the north of the Siwaliks, although this may have been different in the past.

Methods

The procedure followed to collect matched samples of groundwater and aquifer particles has been described [START_REF] Van Geen | Testing groundwater for arsenic in Bangladesh before installing a well[END_REF][START_REF] Van Geen | A transect of groundwater and sediment properties in Araihazar, Bangladesh: Further evidence of decoupling between As and Fe mobilization[END_REF]. Drillers in the Terai region of Nepal use the same drilling method as in Bangladesh and West Bengal. In Vietnam, however, water is forced into sections of PVC pipe with two electric pumps to flush sediment out of the hole instead of using the hand as a check-valve to rhythmically pump sediment-laden water out through the pipe. Drillers in Vietnam also use sections of thin-walled PVC that are cut with a hacksaw, softened over a flame, and glued together at every deployment rather than using threaded metal connectors, which is common in Bangladesh. The central opening of the metal cutter used in Vietnam was modified to allow the needle-sampler to pass through it.

As in the past, the small headspace remaining above the groundwater and sediment slurry in each needle-sampler tube was flushed with N 2 immediately after collection in Vietnam and Nepal. After allowing the sediment to settle for a few minutes, 5-10 mL of groundwater was gently pressure-filtered under N 2 through 0.4 μm syringe-filters and into acid-cleaned scintillation vials with PolySeal caps.

Several days before analysis for As, Fe, Mn, S, Na, Mg, Ca and K by High Resolution Inductively Coupled Plasma-Mass Spectrometry (HR ICP-MS) at LDEO, filtered groundwater was acidified to 1% Optima HCl in the laboratory [START_REF] Cheng | Rapid multi-element analysis of groundwater by highresolution inductively coupled plasma mass spectrometry[END_REF]. The delay in acidification, which reduces the chance of contamination and eliminates the need for shipping strong acid, has been shown to entirely re-dissolve any precipitates that could have formed [START_REF] Van Geen | Monitoring 51 deep community wells in Araihazar, Bangladesh, for up to 5 years: Implications for arsenic mitigation[END_REF].

On the evening of each sampling day, sediment kept in the dark and under groundwater after flushing the headspace of the needle-sampler tube with N 2 was subjected to two different extractions: (1) a hot 10% HCl leach for 30 min [START_REF] Horneman | Decoupling of As and Fe release to Bangladesh groundwater under reducing conditions. Part I: Evidence from sediment profiles[END_REF] and ( 2) an extraction at room temperature in a N 2 -purged 1 M NaH 2 PO 4 solution (pH~5) for 24 h [START_REF] Keon | Validation of an arsenic sequential extraction method for evaluating mobility in sediments[END_REF][START_REF] Jung | Enhanced recovery of arsenite sorbed onto synthetic oxides by l-ascorbic acid addition to phosphate solution: calibrating a sequential leaching method for the speciation analysis of arsenic in natural samples[END_REF]. The acid leachate was analyzed that same evening for dissolved Fe(II) and Fe(III) by colorimetry using ferrozine and later for As by HR ICP-MS. The phosphate extracts were passed through a 0.4 μm syringe filter after 24 h and also analyzed later for As by HR ICP-MS.

Results

This contribution focuses on the geochemical relationships between groundwater and aquifer particles at new sampling sites in Vietnam and Nepal in the context of previous observations from Bangladesh and India based on the same sampling methods. Site-specific spatial patterns in aquifer composition and how these may be controlled by geological history and hydrology is discussed elsewhere in the case of Vietnam [START_REF] Berg | Hydrological and sedimentary controls leading to arsenic contamination of groundwater in the Hanoi Area, Vietnam: The impact of iron-arsenic ratios, peat, river bank deposits, and excessive groundwater abstraction[END_REF]Eiche et al., 2008) and are under investigation in Nepal.

Concentrations of As in filtered groundwater from grey aquifers of Van Phuc village in Vietnam range widely (200±130 μg/L; n=25), with several samples at one end of the spectrum containing ~10 μg/L and several samples approaching 500 μg/L at the other extreme. Concentrations of P-extractable As (2.3±1.5 mg/kg; n=24) and HCl-extractable As (2.5±2.3 mg/kg; n=23) are comparable and vary in concert at this location, but are not related in any systematic way to dissolved As (Fig. 2a,b). Concentrations of dissolved Fe in Van Phuc span more than two orders of magnitude (<0.1 to 40 mg/L), with most groundwater samples containing >100 μg/ L As associated with >3 mg/L dissolved Fe (Fig. 2c). An advanced state of reduction in the particulate phase as well as groundwater anoxia are indicated by leachable Fe(II)/Fe ratios >0.7 and dissolved Mn>0.2 mg/L, respectively, in all sampled intervals (Figs 2c,e). Comparing dissolved S concentrations in groundwater with typical SO 4 2-levels in the Red River (~5 mg/ L, [START_REF] Berg | Hydrological and sedimentary controls leading to arsenic contamination of groundwater in the Hanoi Area, Vietnam: The impact of iron-arsenic ratios, peat, river bank deposits, and excessive groundwater abstraction[END_REF], a likely source of recharge for these shallow aquifers, suggests considerable SO 4

2-reduction has taken place in at least half the groundwater samples from Van Phuc (Fig. 2e).

Whereas higher levels were expected based on a previous survey of Ramgram in Nepal, dissolved As concentration turned out to be relatively low (21±20 μg/L; n=12), with only a single groundwater sample collected with the needle-sampler slightly exceeding 50 μg/L (Fig. 2b). Excluding a (different) interval containing 15 mg As/kg, average P-extractable As concentrations in the particulate phase (0.9±0.8 mg/kg; n=11) are only a factor of two lower than in Van Phuc. There is no clear relation between As concentrations in the solid phase obtained by the two types of sediment extractions, and the P-extractable As content from 3 intervals actually exceeds more than 3 times the level of HCl-extractable As (Fig. 2a). A more systematic distinction between Vietnam and Nepal is that Fe concentrations exceeded 1 mg/L in only one sample from Ramgram (Fig. 2c). Possibly less strongly reducing conditions are also indicated by leachable Fe(II)/Fe ratios <0.5 in 3 samples from Nepal, even if dissolved Mn concentrations of 0.2-0.7 mg/L indicate that groundwater from these intervals is still anoxic (Figs. 2d,f). As in Van Phuc, dissolved S levels <1 mg/L in more than half the groundwater samples indicates significant SO 4 2-reduction in most shallow aquifers that were sampled in Ramgram (Fig. 2e).

Discussion

The new data for dissolved and P-extractable As concentrations in shallow aquifers of Vietnam and Nepal do not fit the relationship previously documented for shallow aquifers of Bangladesh (Fig. 2b). Applying the same criteria as in Bangladesh to identify insufficiently reducing conditions (i.e. diss. Fe <0.2 mg/L and As <10 μg/L) excludes one sample from Vietnam and 5 from Nepal, but even then dissolved As concentrations at a given level of P-extractable As for a considerable number of samples are up to an order of magnitude lower than previously observed in Bangladesh. The needle-sampler data from India also deviate significantly, even after about half the samples from the Hooghly River floodplain are excluded on the basis of either low Fe(II)/Fe ratios or a combination of low dissolved Fe and low dissolved As (Fig. 2b). Unlike the situation in Bangladesh, a single distribution coefficient does not adequately describe the equilibration between dissolved As and P-extractable As in the particulate phase in India, Vietnam, and Nepal.

It is not clear why concentrations of As in groundwater are roughly proportional to the Pextractable As content of aquifer particles across Bangladesh but not at the more recently sampled locations. One notable feature of the combined set of observations is that the relationship established in Bangladesh appears to place an upper bound on the concentration of dissolved As expected for a given level of P-extractable As (Fig. 2b). One possible explanation is that the 1 M phosphate extraction, which constitutes merely an operational definition of the pool of exchangeable As, may be too aggressive outside Bangladesh. Laboratory experiments have shown that even As adsorbed to Fe(III) oxyhydroxides is remobilized in 1 M phosphate because of the strong affinity of P(V) for the same sites [START_REF] Jung | Enhanced recovery of arsenite sorbed onto synthetic oxides by l-ascorbic acid addition to phosphate solution: calibrating a sequential leaching method for the speciation analysis of arsenic in natural samples[END_REF]. [START_REF] Thoral | XAS study of iron and arsenic speciation during Fe(II) oxidation in the presence of As(III)[END_REF] also showed that the structure of Fe(II) oxyhydroxides is limited to very small clusters in the presence of As(III) during Fe oxidation experiments. This suggests that extraction in 1 M phosphate could also release As bound deeply within the natural Fe oxyhydroxides coatings of aquifer particles. It can only be speculated that a milder extraction procedure might be more appropriate for quantifying readily desorbable As and could perhaps have produced a more consistent relationship across a wider range of settings. An alternative explanation is that dissolved As levels are not controlled by adsorptive exchange with an identifiable pool in the particulate phase, even when considering exclusively Fe and SO 4 reducing environments. Differences in mineralogy and/or affinity of surface sites for As are another possibility, but it is not obvious why the mineralogy of sediment reaching different regions of Bangladesh should be any less heterogeneous than sediment deposited in India, Vietnam and Nepal. Some of the differences in the relative efficiency of HCl and P-extractions documented here are consistent with previous observations and may therefore reflect regional trends worth exploring further. BGS/DPHE (2001) reported that oxalate extractions of aquifer sands, which can be considered to target a fraction equivalent to that mobilized by the 1 M phosphate extractions, release only 10-50% of As in the sediment measured by total digestion. [START_REF] Swartz | Mobility of arsenic in a Bangladesh aquifer: Inferences from geochemical profiles, leaching data, and mineralogical characterization[END_REF] and [START_REF] Zheng | Geochemical and hydrogeological contrasts between shallow and deeper aquifers in two villages of Araihazar, Bangladesh: Implications for deeper aquifers as drinking water sources[END_REF] also noted that 1 M phosphate extractions release <30% of total As present in shallow grey sediment at several locations in Bangladesh. Within the suite of samples collected with the needle-sampler in India and Bangladesh, P-extracted As concentrations in sufficiently reducing intervals are on average 0.4-0.5 times lower than HClextractable As levels, respectively (Table 2). It therefore appears that 1 M phosphate extractions systematically release a modest fraction of HCl-extractable As present in aquifer solids from Bangladesh -and this phase appears to be in adsorptive equilibrium with As present in groundwater. [START_REF] Berg | Hydrological and sedimentary controls leading to arsenic contamination of groundwater in the Hanoi Area, Vietnam: The impact of iron-arsenic ratios, peat, river bank deposits, and excessive groundwater abstraction[END_REF], in contrast, report that a similar P-based extraction releases >50% of As in the sediment at several locations in Vietnam. Phosphate extractions also leaves little residual As in the solid phase extractable in HCl for aquifer solids collected with the needle-sampler in Vietnam, as suggested by the comparable P-and HCl-extractable As levels (Table 2). The average P-to HCl-extractable As ratio for aquifer solids from Nepal is similar to the one measured in Vietnam, but more variable and less well constrained. For reasons that are presently not understood, P-extractable As concentrations are significantly lower than HClextractable As levels in two regions that are most affected by elevated As groundwater levels, Bangladesh and West-Bengal, India, whereas the two extractions release comparable amounts of As from the solid phase in two regions where the problem appears to be less widespread. It is speculated that the underlying cause may be related to the extent of transformation of Ascontaining phases in aquifer particles. Observations from Vietnam and Nepal may indicate a relatively mature system, whereas the significant fraction of As contained in the solid phase in the Bangladesh and West Bengal that is not mobilized in 1 M phosphate may indicate incomplete transformation and, consequently, further delay in flushing. One practical implication is that a milder form of extraction, possibly the less aggressive extraction carried out in 0.05 M phosphate by [START_REF] Postma | Arsenic in groundwater of the Red River floodplain, Vietnam: Controlling geochemical processes and reactive transport modeling[END_REF], might be more appropriate for identifying the exchangeable fraction of As in aquifer material across the entire affected area.

Despite the scatter, the new data can be used to place some constraints on the rate at which shallow aquifers are likely to be flushed of their mobilizable As content, assuming it is relatively uniform initially. The simple one-dimensional advection-diffusion model recently presented by van Geen et al. ( 2008) is used with one modification: the distribution coefficient K d is varied between 4 and 40 mL/g to accommodate the new observations (Fig. 2a). Assuming that dissolved and particulate-phase As concentrations are initially equilibrated along the entire flow path at a concentration C i and that recharge water entering the aquifer at one end of the flow path contains no As, the solution of the governing equation is:

where C is the concentration of As in the dissolved phase, x the distance along the flow path, v the advection velocity, D is a hydraulic dispersion coefficient, and R the retardation factor, which is related to the distribution coefficient K d by the expression R = 1 + ρ K d /θ where ρ is the aquifer bulk density, K d is the ratio of the As concentration in the particulate phase on per mass of sediment divided by the As concentration in the dissolved phase, and θ is the aquifer porosity [START_REF] Van Genuchten | Analytical solutions of the one-dimensional convective-dispersive solute transport equation[END_REF].

The focus is on the sensitivity of the model between distances of 50-500 m, i.e. within the range of length of likely flow paths that lead from an area of recharge to the depth range of shallow aquifers and eventually to a discharge area, typically a local stream [START_REF] Stute | Hydrological control of As concentrations in Bangladesh groundwater[END_REF][START_REF] Aziz | Impact of local recharge on arsenic concentrations in shallow aquifers inferred from the electromagnetic conductivity of soils in Araihazar, Bangladesh[END_REF]. Only the central scenario based on a groundwater flow velocity 5 cm/ d (5.8 × 10 -7 m/s) and a dispersion coefficient D = 5.8 × 10 -6 m 2 /s are considered here. Justification for these parameters and model sensitivity to these choices were discussed previously by [START_REF] Van Geen | Flushing history as a hydrogeological control on the regional distribution of arsenic in shallow groundwater of the Bengal Basin[END_REF]. For K d of 4 and 40 mL/g, retardation factors R = 36.8 and 359, respectively, are calculated from an aquifer bulk density of 2.24 g/cm 3 , derived in turn from a porosity of 0.25 and a particle density of 2.65 g/cm 3 .

The model assumes a uniform initial P-extractable As concentration of 4 mg/kg, which corresponds to dissolved As concentrations of 1000 and 100 μg/L for K d = 4 and 40 mL/g, respectively. Under conditions of greatest As mobility, representative of shallow aquifers in Bangladesh, the model shows that a considerable portion of the 50-500 m range of distances along the flow path is already flushed by 0.5 ka and entirely so by 5 ka (Fig. 3). Whereas As concentrations start at a lower level under the alternative scenario associated with a 10-fold higher K d , flushing of the aquifer is also considerably slower. The model predicts elevated As concentrations across the 50-500 distance range, and significant but incomplete flushing of the same shallow aquifers after 5 ka. The scatter in the needle-sampler data suggests that the reality under natural conditions is probably somewhere in between these two extreme scenarios, assuming all sampled deposits are of Holocene age (this has been confirmed in Vietnam and Nepal on the basis of optically stimulated luminescence dating, B. Weinman and A. Singhvi, pers. comm.). The only upside to the fact that As appears to be particularly mobile in groundwater of Bangladesh is that shallow aquifers are likely to be flushed more rapidly in Bangladesh compared to India, Vietnam and Nepal. This is, however, not relevant in terms of human exposure since, hopefully, much of the population in affected regions will no longer consume untreated groundwater within the next few decades. Flushing would be delayed further if, in addition, the fraction of As in the solid phase released in hot HCl is gradually converted to exchangeable As.

No claim is made here that the mechanisms underlying the observed relationships between dissolved and particulate phase As are understood. The new results merely provide some constraints on the range of situations that a more mechanistic model will eventually have to represent. Such a model, possibly based on a more accurate measure of the exchangeable As fraction in the particulate phase, will have to explain why adsorption equilibrium appears to control groundwater As concentrations in Bangladesh more systematically than at the other locations that were sampled. From a practical perspective, it would also be important to understand under what conditions the considerable pool of P-extractable As stored in the particulate phase in less reducing aquifers characterized by leachable Fe(II)/Fe ratios<0.5 or dissolved Fe <0.1 mg/L could potentially be released to groundwater in the future.

Conclusions

The proportionality between dissolved and P-extractable As concentrations documented for shallow aquifers of Bangladesh does not appear to hold in other more recently sampled regions of southern Asia. Even after excluding intervals where conditions may not be sufficiently reducing, As concentrations in groundwater of India, Vietnam and Nepal are frequently up to an order of magnitude lower at given levels of P-extractable As than in Bangladesh. One possible explanation is that the 1 M phosphate extraction releases As from a less labile fraction in the particulate phase that, for reasons presently unclear, is more important outside Bangladesh. To a first approximation, much of the needle-sampler data obtained from grey aquifers throughout the affected region can nevertheless be ascribed to adsorptive equilibration according to a distribution coefficient that ranges from 4 to 40 mL/g. Using this formulation and a simple model of groundwater flow, it is estimated that late Holocene aquifers have been, or will be, flushed of much of their initial As content 0.5 to 5 ka after their deposition. The uncertainty of such estimates is considerable. More robust is the conclusion that the elevated As content of groundwater throughout southern Asia is probably a transient feature, although on a time scale that extends considerably beyond the lives of the presently exposed population. 2008) with additional well As data for Nepal from [START_REF] Shrestha | The State of Arsenic in Nepal-2003, National Arsenic Steering Committee[END_REF]. In the case of Vietnam, the 3 shaded areas correspond to high, medium and low risk of groundwater exceeding 50 μg/L As, as inferred from geological information by [START_REF] Berg | Magnitude of arsenic pollution in the Mekong and Red River Deltas-Cambodia and Vietnam[END_REF]. The inset is an expanded view of the Red River delta area. The shape of the symbols corresponds to the symbols used to distinguish different locations in Fig. 2. The map was created from the Digital Atlas of South Asia of [START_REF] Hearn | Global GIS Database: Digital Atlas of South Asia[END_REF].

Fig. 1 .

 1 Fig. 1. Provenance of needle-sampler material from Bangladesh, India, Vietnam and Nepal shown on a map of As-affected areas in southern Asia. The majority of tubewells sampled and analyzed by Dr. Dipankar Chakraborti's group at Jadavpur University located within red (dark grey in the printed copy) and light blue (lightest grey) areas in India and Bangladesh contain >50 μg/ L and <50 μg/L As, respectively. Orange (intermediate grey) areas indicate a mixed distribution of high-and low-As wells. Modified from Métral et al. (2008) with additional well As data for Nepal from[START_REF] Shrestha | The State of Arsenic in Nepal-2003, National Arsenic Steering Committee[END_REF]. In the case of Vietnam, the 3 shaded areas correspond to high, medium and low risk of groundwater exceeding 50 μg/L As, as inferred from geological information by[START_REF] Berg | Magnitude of arsenic pollution in the Mekong and Red River Deltas-Cambodia and Vietnam[END_REF]. The inset is an expanded view of the Red River delta area. The shape of the symbols corresponds to the symbols used to distinguish different locations in Fig.2. The map was created from the Digital Atlas of South Asia of[START_REF] Hearn | Global GIS Database: Digital Atlas of South Asia[END_REF].

Fig. 2 .

 2 Fig. 2. Scatter plots of groundwater and aquifer particle properties for material collected with the needle-sampler in Bangladesh, India, Vietnam and Nepal. Note the logarithmic scale of all properties, with the exception of leachable Fe(II)/Fe ratios. Open symbols identify samples collected from insufficiently reducing samples, as identified in van Geen et al. (2008) by either leachable Fe(II)/Fe <0.5 or dissolved Fe concentrations <0.2 mg/L in samples with dissolved As <10 μg/L. Solid and dashed lines in (a) correspond to equal and 10-fold higher As concentrations in HCl-extractions, respectively, compared to As levels in P-extractions. Solid line and dashed lines in (b) correspond to K d = 4 and 40 cm 3 /g, respectively. The arrow in (d) points to a single sample from Bangladesh containing an unexplained and extremely elevated level of 190 mg S/L.

Fig. 3 .

 3 Fig. 3. Predicted distribution of groundwater As concentrations along a flowpath according to the advection-diffusion model described in van Geen et al. (2008) after 50, 500 and 5000 a for K d s of 4 and 40 cm 3 /g. Note the logarithmic scale of distance along the flow path. Also shown is the likely 50-500 m range of lengths of flow paths from recharge to discharge areas within shallow aquifers.

Appl Geochem. Author manuscript; available in PMC 2009 November 1.

Acknowledgments

Expansion of arsenic studies beyond Columbia University's initial study area of Araihazar, Bangladesh, was supported by NSF grant EAR 0345688. Support over the long term has been provided by NIEHS Superfund Basic Research Program grant P42 ES10349. The LGIT-OSUG field work in Nepal was supported by the EC2CO-INSU-CNRS program. We are grateful to numerous colleagues, students, and drillers for their contributions to the field work, including Linda Smith of Filters for Families in Nepal. The manuscript benefited from insightful comments provided by George Breit and Flemming Larsen. This is Lamont-Doherty Earth Observatory contribution number XXXX.