Pavel Karas
email: xkaras1@fi.muni.cz

Vincent Morard
email: vincent.morard@mines-paristech.fr

Jan Bartovsky
email: j.bartovsky@esiee.fr

Thierry Grandpierre
email: t.grandpierre@esiee.fr

Eva Dokladalova
email: e.dokladalova@esiee.fr

Petr Matula
email: petr.matula@mines-paristech.fr

Petr Dokládal
email: petr.dokladal@mines-paristech.fr

Jan Bartovský

Eva Dokládalová

GPU Implementation of Linear Morphological Openings with Arbitrary Angle

Keywords: morphology, opening, closing, linear, 1-D SE, parallel, efficient, algorithm, implementation, GPU

HAL is

tistical measures called granulometries [START_REF] Vincent | Fast opening functions and morphological granulometries[END_REF][START_REF] Batman | Heterogeneous morphological granulometries[END_REF][START_REF] Vincent | Granulometries and opening trees[END_REF][START_REF] Urbach | Connected rotation-invariant size-shape granulometries[END_REF] and of a class of non-linear morphological filters called Alternate Sequential Filters (ASF) [START_REF] Serra | An overview of morphological filtering[END_REF][START_REF] Heijmans | A new class of alternating sequential filters[END_REF].

In practical applications, the granulometries allow estimation of a priori unknown geometrical characteristics of objects in the image. For illustration, we can cite (1) medical imaging applications e.g., blood cell classification [START_REF] Theera-Umpon | Counting white blood cells using morphological granulometries[END_REF], (2) automated document analysis [START_REF] Bagdanov | Granulometric analysis of document images[END_REF], or (3) industrial control [START_REF] Outal | A new method for estimating the 3d size-distribution curve of fragmented rocks out of 2d images[END_REF]. The role of the ASF filters is to reduce the noise while preserving the principal features in the image. They represent the principal element of numerous applications e.g., texture analysis [START_REF] Talukdar | Estimation of fractal dimension using alternating sequential filters[END_REF] or remote sensing [START_REF] Sigurjonsson | Street tracking based on sar data from urban areas[END_REF]. Even the morphological openings themselves are useful for their filtering properties in some industrial applications such as [START_REF] Kowalczyk | Application of mathematical morphology operations for simplification and improvement of correlation of images in close-range photogrammetry[END_REF].

Generally speaking, to obtain the desired result-size distribution or filtering effect-we have to use a sequence of openings and/or closings with varying parameters of the applied computing window, so-called structuring element (SE). For a given shape of the SE, these variable parameters are the progressively increasing size of SE and rotation angle. In order to ensure the exhaustivity of the result, applications often require computing of an enormous number of iterations with greater SE, often approaching hundreds of pixels. Considering continually increasing image resolution used in industrial applications, one can intuitively feel that it results in overwhelming requirements on the computing power. This is true even despite recent efficient algorithms [START_REF] Morard | Linear openings in arbitrary orientation in O(1) per pixel[END_REF][START_REF] Bartovský | Fast streaming algorithm for 1-d morphological opening and closing on 2-d support[END_REF].

In this context, we study how to efficiently implement the above mentioned operators on graphics cards with the objective to reduce these computing requirements on the system. Initially, graphics cards were designed for graphics purposes only and were not programmable. Based on numerous parallel processors they were very powerful compared to their price. Current GPUs have passed the one Tera FLOPS barrier, and there is no need to use dedicated graphics languages any more since several frameworks have been developed for GPGPU 1 purposes: CUDA [START_REF]NVIDIA GPU Computing Developer Home Page[END_REF] by nVidia and OpenCL [START_REF] Ouzounis | OpenCL[END_REF] by the Khronos Group are today the most popular in the GPGPU community. Both are based on C language extensions. Notice that in this paper, we use the CUDA language for all presented implementations on GPU. Nevertheless, the principles remain the same when moving to a different language.

In our study, we focus on the morphological openings using linear SE under arbitrary angle, essential in a wide range of practical applications dealing with linear image structures such as fingerprint analysis, geosciences [START_REF] Obara | Identification of transcrystalline microcracks observed in microscope images of a dolomite structure using image analysis methods based on linear structuring element processing[END_REF], and vessel segmentation [START_REF] Zana | Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation[END_REF].

The main novelty of this paper consists of an efficient and optimized GPU implementation of the algorithm by Bartovsky et al. [START_REF] Bartovský | Fast streaming algorithm for 1-d morphological opening and closing on 2-d support[END_REF], which turned out to be the most suitable for a parallel implementation on GPU. We also present performance comparisons and experimental results of all implementations on real data. It turned out that the proposed GPU implementation can compute linear openings for 180 directions of an image of size 640×640 pixels within 60 ms.We also show that the proposed implementation is significantly faster than the state-of-the-art implementation in the OpenCV_GPU library [START_REF] Garage | OpenCV_GPU[END_REF], especially for larger SEs.

In the remainder of this paper we start by the introduction of the mathematical background in Section 2. Section 3 reviews and compares the state of the art of efficient implementations for computing linear morphological openings of different orientations. Section 4 presents the first GPU implementation of two candidate algorithms by Bartovsky [START_REF] Bartovský | Fast streaming algorithm for 1-d morphological opening and closing on 2-d support[END_REF] and Morard [START_REF] Morard | Linear openings in arbitrary orientation in O(1) per pixel[END_REF]. Afterwards, the Section 5 describes the optimized implementation of Bartovsky algorithm, including the parallelism enhancement discussion. Section 6 illustrates the use of this efficient implementation in practical applications. It includes also the discussion of overall experimental results. Finally, the conclusions recall the main contributions of our work and briefly introduces its perspectives.

Basic Notions

Before defining morphological openings and closings, we define erosions and dilations for gray-scale images. Let a mapping f : Z 2 → R be a gray-scale discrete image, and F be the collection of all such images. The support used throughout this paper shall be a rectangular subset domain of Z 2 .

The erosion and dilation are mappings ε B , δ B : F → F , parameterized by the so-called structuring element (SE) B, B ⊂ Z 2 . Here, we restrict the family of SE to flat and translation-invariant. Hence, the dilation and erosion are, as usually, defined by

δ B ( f ) = b∈B f b ; ε B ( f ) = b∈ B f b (1)
where the hat denotes the transposition of the structuring element, i.e., B = {x | -x ∈ B}, and f b denotes the translation of the function f by some vector b, and and denote the supremum and infimum on a collection of functions. Similarly, the morphological openings and closings are mappings γ B , ϕ B : F → F , also parameterized by B. The transposition B used in the definition of erosion ensures that the dilation and erosion form a so-called adjunction pair. This allows us to obtain the morphological opening and closing by concatenation of the dilation and erosion:

γ B ( f ) = δ B [ε B ( f )] ; ϕ B ( f ) = ε B [δ B ( f )]. (2)
Erosions and dilations are dual under complementation, i.e. for functions δ B = -ε B (f ), see e.g. [START_REF] Soille | Morphological Image Analysis: Principles and Applications[END_REF]. Consequently, γ B and ϕ B are also dual and all algorithms developed for openings can easily be used also for closings.

In the sequel, we focus on linear morphological openings and closings obtained with SE in a form of a 1-D digital line segment of the length λ and orientation α, denoted by γ α λ and ϕ α λ , respectively. Further operators can be built based on linear openings. The first operator, ∨γ λ ( f ), is computed by taking the supremum of the openings by digital line segments in all orientations; the second operator, ζ λ ( f ), can extract the local orientation of linear structures:

∨γ λ ( f ) = α∈[0,180] γ α λ ( f ) ; ζ λ ( f ) = arg α∈[0,180] γ α λ ( f ). (3)
As an example, the effect of these operators on a real image is presented in Fig. 1.

Opening Algorithms

Opening algorithms can be divided into three classes: erosion and dilation chaining (a two-stage algorithm), direct computation, and algorithms based on connected component tree building. Two-stage algorithm: computed by using Eq. ( 2), it is the simplest and historically earlier approach used to compute openings. Noting N the number of pixels of the image and L the size of the SE, the complexity of the naïve approach is O(N × L). The computational complexity was improved throughout the years to make the algorithm feasible for practical applications. In 1985, Pecht [START_REF] Pecht | Speeding-up successive minkowski operations with bit-plane computers[END_REF] defined a logarithmic decomposition of the SE. This decomposition, which removes most of the redundancy, was further extended by Coltuc [START_REF] Coltuc | On fast running max-min filtering[END_REF], reducing the complexity to O(N × log L). Later, the complexity was further reduced to linear O(N), hence independent of the size of the structuring element, by Van Herk, and Gil and Werman [START_REF] Van Herk | A fast algorithm for local minimum and maximum filters on rectangular and octagonal kernels[END_REF][START_REF] Gil | Computing 2-d min, median, and max filters[END_REF]. The linear algorithm will be referred to as the HGW algorithm hereafter. In Direct Computation: for the family of openings with 1-D SE, new algorithms were introduced recently to directly compute openings in only one scan of the entire image, preserving the complexity of O(N). In [START_REF] Van Droogenbroeck | Morphological erosions and openings: fast algorithms based on anchors[END_REF], Van Droogenbroeck and Buckley introduced an algorithm based on anchors, allowing very fast computation of the linear openings. The anchors are points, which are not affected by the opening. Nevertheless, the algorithm uses a histogram. The main drawback of using a histogram is the increasing memory consumption for finely quantized data. Even though the memory is no longer an issue for PC architectures, it becomes a penalizing factor for parallelized implementations on other architectures with limited memory like GPU. Secondly, search over a long histogram becomes costly and finely, for floating point accuracy, the histogram can not be used at all. Later, Morard et al. [START_REF] Morard | Linear openings in arbitrary orientation in O(1) per pixel[END_REF] introduced a very simple algorithm based on an ordered stack of cords where a cord refers to a continuous set of pixels of intensity greater than or equal to a certain value I. With the inclusion relation between cords and by analyzing the length of each cord, the computation of linear openings and of linear granulometries is straightforward. Finally, Bartovsky et al. [START_REF] Bartovský | Fast streaming algorithm for 1-d morphological opening and closing on 2-d support[END_REF] also developed an algorithm to build linear openings. It sequentially scans the signal and erases every peak narrower than SE. Further explanations on these algorithms are given in section 4.

Connected component tree: such approach was described in [START_REF] Garrido | Extensive operators in partition lattices for image sequence analysis[END_REF]. The approach is based on building connected components, hence it can be adopted for more complex tasks such as watershed segmentation [START_REF] Matas | Parallel algorithm for concurrent computation of connected component tree[END_REF]. The drawback of the algorithm is that the complexity depends on the number of gray levels in the image and requires random access data, consequently it is not adapted for applications running under strict time constraints.

Parallel Implementations

There are several implementations of HGW algorithm in the literature since it can be easily parallelized. Brambor [START_REF] Brambor | Algorithmes de la Morphologie Mathématique pour les architectures orietées flux[END_REF] described a parallel implementation of the HGW algorithm on SIMD architectures. Their implementation was tested on an Intel CPU with the SIMD-SSE2 instruction set. Clienti [START_REF] Clienti | An efficient hardware architecture without line memories for morphological image processing[END_REF] improved the HGW algorithm and implemented2 it on a SIMD architecture as well. Domanski et al. [START_REF] Domanski | Parallel van Herk/Gil-Werman image morphology on GPUs using CUDA[END_REF] used CUDA to implement the HGW algorithm on GPU, achieving 13-33 × speedup. There are several drawbacks of the algorithm as it computes openings and closings by dilation-erosion chaining, which requires more computations than direct approaches. It also has larger memory requirements [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF].

On the contrary, there are few parallel implementations of the component tree algorithms, among which we can cite Wilkinson [START_REF] Wilkinson | Concurrent computation of attribute filters on shared memory parallel machines[END_REF], Menotti-Gomes [START_REF] Menotti-Gomes | 1D Component tree in linear time and space and its application to gray-level image multithresholding[END_REF] on multicore, and also Matas [START_REF] Matas | Parallel algorithm for concurrent computation of connected component tree[END_REF] on ccNUMA 4-core. They are effectively so complex that it is difficult to exhibit some parallelism in these complex algorithms.

Finally, in order to improve the computing efficiency by parallel implementation, direct computation algorithms seem to be the best candidates compared to HGW and component-tree algorithms. This is why we focus on this class of algorithms in the remainder of this paper.

Selection of a Direct Linear Opening Algorithm

We need to select the best sequential algorithm candidate for a parallel implementation leading to the most efficient execution on a GPU platform. As explained above, such an algorithm has to allow arbitrary angles computing using direct linear opening for lower computation and complexity requirements. Hence, there are three algorithms available to benchmark and compare: (1) Van Droogenbroeck [START_REF] Van Droogenbroeck | Morphological erosions and openings: fast algorithms based on anchors[END_REF] referred to as VanD, (2) Morard et al. [START_REF] Morard | Linear openings in arbitrary orientation in O(1) per pixel[END_REF] referred to as Morard and (3) Bartovsky et al. [START_REF] Bartovský | Fast streaming algorithm for 1-d morphological opening and closing on 2-d support[END_REF] referred to as Bartovsky. To this list, we will add Clienti algorithm [START_REF] Clienti | An efficient hardware architecture without line memories for morphological image processing[END_REF] although it is not a direct computation based algorithm. Effectively, this is one of the fastest HGW implementation that can consequently give a useful comparative point. It will be referred to as Clienti.

For a fair comparison, all these algorithms were implemented using the same image-processing library with the same interface and with all the optimization flags turned on. We used only one core of an Intel Core i7-870 2.93 GHz CPU for this benchmark. These algorithms have been applied on the texture images shown in Fig. 12(a),(b).

Execution performances of the four algorithms for both horizontal and vertical openings are presented in Fig. 2. The performance P is computed as P = N/t, where N is the image size and t is the computation time. This measure is consequently independent of the image size. Note, however, that the performance actually does depend on the image size (see e.g. Fig. 2(b)). Each marked value in the plots represents the performance computed from the mean computation of 100 openings. For each image the length of a structuring element was chosen to be equal to 5 % of the image width/height because, in most applications, the length of SE is considerably smaller than the image size.

From the benchmarks we see that for vertical SEs and for large images, the performances significantly decreases for all algorithms, except for the Bartovsky. This is explained by the fact that this algorithm accesses the data sequentially even for vertical structuring elements. It is clear that VanD algorithm is the fastest, followed by Morard, Bartovsky, and Clienti, for both vertical and horizontal orientation. While VanD achieves the best performance, it is nevertheless unsuitable for parallelization on a GPU because of its high memory requirements especially for higher bit depths and arbitrary oriented SEs. In contrast, Morard and Bartovsky algorithms are able to compute openings and closings efficiently regardless the orientation of the SE or the data type precision with minimum memory requirements. Therefore, we selected these two algorithms for parallelization and implementation onto the targeted GPU architecture.

Algorithm 1: Morard algorithm: G ← Open1D(F, λ , S) Input: F -input 1-D signal, λ -size of SE, S -

Basic Implementation on GPU

Morard and Bartovsky Algorithms

In this section, we briefly introduce the Morard and Bartovsky algorithms for computation of morphological openings and closings with a linear SE. The detailed description can be found in [START_REF] Morard | Linear openings in arbitrary orientation in O(1) per pixel[END_REF][START_REF] Bartovský | Fast streaming algorithm for 1-d morphological opening and closing on 2-d support[END_REF]. Both algorithms are able to compute the operation in O(N) time with respect to the image size and O(1) with respect to the size of SE. The Bartovsky algorithm was originally designed for streaming architectures such as FPGA and hence performs scanning of the input data in sequential order. Nevertheless, it can be simply modified to perform scanning along lines according to the orientation of SE, much like the Morard algorithm. During the scan, intensity and position of each pixel is stored in an auxiliary data structure if necessary. Whereas the Morard algorithm uses the LIFO stack, the Bartovsky algorithm uses the FIFO queue. For arbitrary directions, both algorithms use the Bresenham's lines, as described in [START_REF] Soille | Recursive implementation of erosions and dilations along discrete lines at arbitrary angles[END_REF].

In the Morard algorithm, the image line is scanned for pixels where the intensity changes. Whenever the current pixel intensity is higher than the preceding, the current pixel is pushed to the stack. In the opposite case, the stack is being emptied while necessary. The algorithm needs to store an extra bit for a boolean flag indicating the status of a pixel. The size of the stack is limited only by the size and the bit depth of the image. After processing the stack, the output is written. The outputs are irregular. A pseudo-code is shown in Alg. 1.

In the Bartovsky algorithm, the image line is scanned for so-called peaks. According to a peak configuration, either a pixel is pushed to the queue or the queue is being emptied. The size of the queue is limited by the size of SE. After Algorithm 2: Bartovsky algorithm: G ← Open1D(F, λ , Q)

Input: F -input 1-D signal, λ -size of SE, Q -pointer to FIFO queue Result: G -output 1-D signal Data: Q -a double-end queue of pairs (value, position) initialize Q; for rp ← 0 to |F| -1 do while F[rp] ≤ Q.back(value) do process Q; Q.push({F[rp], rp}); if rp = Q.front(position) + λ then Q.pop(); if rp ≥ λ then G[rp -λ ] ← Q.front(value);
processing the queue, the output is written. The writes are regular and follow the reads with the well-defined latency corresponding to the size of SE. A pseudo-code is shown in Alg. 2.

GPU Parallel Architecture

Before we describe basic GPU implementations of the algorithms, we briefly review some properties of the GPU architecture, as shown in Fig. 3.

From the hardware point of view, a GPU consists of hundreds of so-called streaming multiprocessors (SM). Each includes a number of shader processor (SP) cores, a number of registers, a small shared memory providing communication between SPs, and fast ALU units for hardware acceleration of transcendental functions. One global memory is shared by all SMs and provides a capacity in order of GB and the memory bandwidth in order of 100 GB/s. There are also two additional read-only cached memory spaces accessible by all threads: the constant and texture memory spaces. They can help programmers to improve the performance of their implementations [START_REF]CUDA Toolkit Reference Manual. 2701 San Tomas Expressway[END_REF][START_REF] Kong | Accelerating MATLAB Image Processing Toolbox functions on GPUs[END_REF]. In some recent GPU architectures, such as FERMI by nVidia [START_REF]FERMI Tuning Guide[END_REF], the global memory is cached as well.

From the programmer's point of view, every program consists of two parts, a host code for CPU, and a kernel code for GPU. Before executing a kernel, the host program allocates a memory on GPU and transfers data if necessary. Then a kernel is configured and executed. The configuration defines the number of threads allocated for kernel execution. Generally, the number of threads should be much higher than number of processing units of GPU. This approach allows (1) the proper scaling on various hardware configurations, and (2) hiding the memory latencies. The threads form groups called blocks (as in CUDA [START_REF]CUDA Toolkit Reference Manual. 2701 San Tomas Expressway[END_REF]) or work-groups (as in OpenCL [START_REF]OpenCL Course: Introduction to OpenCL Programming[END_REF]), following the hierarchy of SMs and SPs. Synchronization of threads and data sharing is possible within the group only. The threads are executed concurrently in warps, usually of 32 threads each.

GPU Implementation of Morard and Bartovsky Algorithms

Regarding the design of both algorithms, the input image is scanned in the sequential manner. However, all lines of the image can be scanned concurrently, as shown in Fig. 4. Thus, the parallelism can be introduced by binding individual threads to individual lines of the image. Each image line is processed independently using its own algorithm-dependent stack or queue, respectively. This requires the GPU memory to be large enough to contain all auxiliary data structures. Here, the Bartovsky algorithm is favorable, since the sizes of the queues used are limited to the size of SE whereas the stacks used by Morard algorithm are generally limited by the size of the image.

The mapping of threads for all SE orientations is described in Fig. 4. For an arbitrary angle α, the overall number of threads can be simply computed as follows:

T = w + ⌈h| cot α|⌉, α ∈ [45 • , 135 • ) ∪ [225 • , 315 • ), (4a)
T = h + ⌈w| tan α|⌉, α ∈ [-45 • , 45 • ) ∪ [135 • , 225 • ), (4b)
where w and h are the width and the height of the input image, respectively. Thus, generally T > w or T > h, respectively. Some of threads spend a part of the computation time outside the image domain where they do not compute anything. Thanks to thread locality, this brings little overhead only because in the GPU thread scheduler, thread warps that fall outside the image domain are quickly replaced by those that fall inside.

It should be noted that the processing of both the stack and the queue is data-dependent, thus data accesses are irregular. Therefore, during this part of the computation, threads are divergent. This limits the overall performance of the parallel implementation. In the Bartovsky algorithm, the data accesses to both input and output images are regular.

Experimental Results

The basic GPU implementations, described in the previous section, were compared on nVidia Tesla C2050 GPU with 14 MPs at 1.15 GHz and 3 GB RAM. Again, all experiments were made with the texture images shown in Fig. 12(a), (b). Results are shown in Fig. 5. The presented performance does not include the times needed for the data transfer between CPU and GPU. Contrarily to the performance on CPU, the Bartovsky algorithm achieves better performance than the Morard algorithm. However, its performance is highly sensitive to the SE orientation, as shown in Fig. 5. For a horizontal SE, the performance is significantly lower due to extensive L1 cache misses, as detected by Nsight Visual Profiler [START_REF]NVIDIA Parallel Nsight[END_REF]. As the L1 cache is stored in memory banks and memory accesses are synchronized, most of cache queries are directed to the same memory bank leaving other memory banks unused. This is not the case of the Morard algorithm since it generally keeps the memory accesses irregular.

Efficient Implementation of the Bartovsky Algorithm

The choice of the algorithm to optimize is guided by the analysis of the Bartovsky and Morard algorithms in sections 4.3 and 4.4. Despite higher performances of the Morard algorithm on CPU, the Bartovsky algorithm has several advantages:

(1) Both the data accesses to the input and output image are regular. This helps to make the thread execution synchronous and reduces the thread divergence.

(2) The maximum length of the FIFO queue is limited by the length of SE. This strongly limits the memory requirements. Recall that we are to bind one thread per image line (Fig. 4). For large images, with potentionally many threads, it is important to have a small memory footprint of each thread.

In the basic implementation on GPU (contrarily to CPU), the Bartovsky algorithm is sensitive to the SE orientation. Hence, we shall introduce several optimization steps not only to increase the overall performance but, particularly, to keep the performance stable for all orientations. These steps are described in the following sections.

To prove the choice of the Bartovsky algorithm, we applied the optimization steps also on the Morard algorithm and performed tests to compare both optimized GPU implementations. The results are shown in Section 6.

Parallelism Enhancement

In the basic implementation, the parallelism was introduced by mapping GPU threads to individual image rows or columns, creating the grid of h or w threads, respectively (where h and w are height and width of a 2-D input image, respectively). However, if an input image is not large enough, the GPU's MPs are not fully occupied. Therefore, we introduce more parallelism by splitting the image into two or more parts. In the following text, we refer to them as zones. As each pixel can be affected by (2λ -2) neighboring pixels (where λ is the length of SE), the zones need to overlap by 2(λ -1) pixels. An example of the image split for a vertical SE is shown in Fig. 6.

It is evident that the theoretical speed-up s that can be achieved depends on the SE length λ , the image size and the number of zones Z. For vertical SE we get the following:

s = Z 1 - (Z -1)(2λ -1) h , (5)
where h is the image height. For small λ we get s ≈ Z, while for λ ≈ h/(2Z) we get s ≈ 1. It should be noted that for large input images it is not necessary to introduce more parallelism by splitting the image, it can even decrease the performance. The optimal choice of the parameter Z is discussed in section 5.4.

Optimization of Data Accesses

To reduce the memory latency, the usage of the global memory should be avoided where possible. Alternative memory spaces can be used for both the input image and the FIFO queues. The input image can beneficially be stored in the read-only cached texture memory. This is true also for the recent FERMI architecture, which introduced L1 cache [START_REF]FERMI Tuning Guide[END_REF][START_REF] Nickolls | The GPU Computing Era[END_REF], because the texture memory helps significantly improve the performance for horizontal SEs where L1 cache fails due to bank conflicts. The FIFOs can be stored in the shared memory. As the amount of the available shared memory is limited, the maximum length l max of the FIFO is limited to l max = S/(S b × d), where S is the amount of the shared memory available, S b is the number of threads per block (workgroup) sharing the memory, and d is the size of the data type. Hence, the most recent values of the FIFO are stored in the shared memory in a circular buffer so the position of the first element varies during the computation. The rest of the FIFO is stored in the global memory. The FIFO storage is shown in Fig. 7.

The constraint on the FIFO length mentioned above can be considered as a hard limit. However, using as much shared memory as possible can lead to a performance decrease because in that case, the shared memory is allocated for one thread block per MP only. Since the MPs are capable of execution of more thread blocks, the optimal FIFO length needs to be chosen carefully, as discussed further.

To conclude our choice of memory spaces, using L1 and L2 caches for input data turned out to be inefficient for some SE orientations, as explained in Section 4.4. The texture memory is read-only, cached, and can be allocated up to the size of the device memory. Therefore, it is suitable for the input image. The size of the texture cache is relatively small (8 kB per multiprocessor) but sufficient, thanks to the spatial locality. The shared memory is rewritable, relatively small (up to 48 kB per MP), and fast. In terms of latency, it is comparable with registers, provided that block conflicts are avoided-which is guaranteed in our implementation. Thus, it is suitable for FIFO caching.

Consideration of Arbitrary Orientations

As noted in section 4.4, the performance of Bartovsky algorithm on GPU is sensitive to SE orientation. Thus, a careful attention should be paid to this issue when computing openings in arbitrary directions. By using the texture memory, the input image is cached and the latency of memory reads is reduced. The final step is to optimize the memory writes to the output image. Experiments proved that openings with SEs of orientation α ∈ (-45 • , 45 • ) are computed faster when writing the output image transposed. The output image is subsequently corrected using the modified transpose kernel from [START_REF]CUDA SDK Code Samples[END_REF].

Configuration of Performance Parameters

In the optimized GPU implementation, we have introduced several parameters that influence the performance: the block size (work-group size) S b , the number of zones Z, and the optimal length l of the FIFO buffer allocated in the shared memory. They all depend on these properties: the image dimensions (w or h), the SE length λ , the number of the MPs N MP available on the used device, the shared memory size per MP S, the warp size W and the number of blocks (workgroups) N b that can be executed on a single MP. For optimal configuration, the following set of rules should be satisfied:

1. S b should be multiple of W , 2. Z ≥ (S b × N b × N MP )/T where T is defined in Eq. ( 4), 3. Z ≤ h/(2λ ) or Z ≤ w/(2λ ) following Eq. ( 5), 4. l ≤ S/(S b × d × N b ), 5. l should be a power of two allowing the bitwise operator "&" to be used instead of the costly modulo operator for addressing the FIFO queue items.

It should be noted that some of the rules cannot be satisfied in some cases. In particular, rules ( 2) and ( 3) can be conflicting for small input images. Performances for some parameter configurations are presented in the following section. CUDA programmers are advised to use a tool called "CUDA Occupancy Calculator" which can help to compute the optimal kernel configuration [START_REF]CUDA SDK Code Samples[END_REF].

Experimental Results

We made the performance analysis of the optimized GPU implementations, based on images with clear linear structures (see Fig. 12), and the results of our comparisons are shown in Fig. 8, 9, and 10.

In the first experiment (Fig. 8), the comparison of optimized GPU implementations of Bartovsky and Morard algorithms is shown, in analogy to Fig. 5. We assumed that the former is more suitable for the GPU architecture than the latter. This assumption was confirmed by the experiments. Thus, in the following experiments, we used the more successful implementation.

In the second experiment (Fig. 9 and 10), we compared the performance of our GPU implementation, referred to as "GPU (Bartovsky)", to the corresponding CPU implementation, referred to as "CPU (Bartovsky)", and also to the stateof-the-art implementation in the OpenCV library with the GPU support (so-called OpenCV_GPU) [START_REF] Garage | OpenCV_GPU[END_REF], referred to as "GPU (OpenCV)". It turns out that our GPU implementation is approximately 10-50 × faster than the CPU implementation, depending on the input data size and the length of the SE. Despite the fact that the Bartovsky algorithm itself is sequential so the parallelism introduced in its GPU implementation is limited, our implementation is faster than the OpenCV_GPU in every case. Whereas for small SEs the difference between the two GPU implementations is negligible, for larger horizontal and vertical SEs the speedup is up to 50 ×. For diagonal (and arbitrarily oriented) SEs, the performance of the OpenCV_GPU implementation falls down very quickly. This is because this implementation uses the NVIDIA Parallel Primitives (NPP) library [START_REF]NVIDIA Performance Primitives (NPP) Library User Guide[END_REF] which supports only simple SE shapes, therefore the line SE has to be represented by a corresponding 2-D rectangle, i.e. a matrix with elements correctly set to 0 or 1.

The most important performance limit of our implementation is the number of threads that can be executed. If the Table 1 Optimal kernel configuration for nVidia graphics cards and for a vertical SE. Choice of the optimal block size S b , the Z parameter, and the size of the shared memory l depend on many parameters, as described in section 5.4. The Z parameter should be decreased if the SE size is too large. Compute capability refers to a core architecture of nVidia graphics cards [START_REF]CUDA C Programming Guide[END_REF]. input image is too small, there is not enough threads and the GPU's is underused. This can be avoided by splitting the image in zones (refer to Fig. 6). Adjacent zones need to overlap to avoid border effects. For large SE sizes, the overlap becomes large, with the consequence that the performances decrease, see Fig. 10(a). For larger image sizes, the decrease is proportionally lesser, see Fig. 10(b). Thanks to the optimizations described in section 5.3 we achieved stable performance for all SE orientations. The performance was tested for all α ∈ [0 • , 180 • ), although for the sake of simplicity, the graphs show only three orientations. To conclude, our GPU implementation can be used for an arbitrary SE length and direction, achieving the performance more than 1000 Mpx/s. This allows computing one opening on a 40 Mpx image with any 1-D SE in any orientation.

The contribution of splitting the image into zones is shown in Fig. 11 and compared with a theoretical speedup, as computed by Eq. ( 5). For smaller images, the optimization increases the performance as expected. Note that for large SE sizes, the performance does not decrease as quickly as theoretically predicted. It is because the further parallelism introduced by the split not only occupies more MPs but it also helps to hide memory latencies. For larger images, the number of threads is large enough to occupy MPs, hence the image split does not introduce a further speedup.

The optimal choice of parameters for nVidia graphics cards and for a vertical SE is shown in Table 1. For AMD Radeon graphics cards, the optimum values may differ. For other orientations, the optimal parameters are analogous.

The performance values do not include the data transfers between CPU and GPU. According to our measurements, the time needed to transfer data is comparable with the time needed for the computation of a single opening, hence the overall speedup is half. Here, the GPU implementation is favorable for images larger than 6 Mpx. In the computation of multiple openings, the data transfer overhead is negligible. The performance values for our GPU implementation were obtained on top-class Tesla C2050 GPU (current price 1500 EUR), but we did several test also on 10 × cheaper GeForce GTX 470 GPU, and the results were comparable.

Practical Applications

In practice, linear openings and closings can be used for the detection of either local or global orientations of linear structures. Hence, we tested and compared two CPU implementations (Bartovsky, Morard) and our GPU implementation of linear openings and closings based on images from three different application domains, namely fingerprint analysis, texture characterization, and document analysis. In all cases, we computed a set of linear openings allowing massive parallelism on GPU simply by computing linear openings in all directions concurrently avoiding the problem with insufficient MPs occupancy. We will show that this leads to a significant speed-up even for small input images. The benchmark results for all applications are shown in Table 2.

Local Orientation of Linear Structures in Fingerprint Images

The first application was the computation of local orientation of linear structures in a fingerprint image (Fig. 1). The operator ζ λ ( f ) was computed according to Eq. ( 3). The GPU implementation achieves a significant speedup (approximately 27 ×) even for small input images (0.9 megapixels).

Angular Spectrum of Texture Images

The second application was the computation of the angular spectrum of texture images. The spectrum was computed in order to find the most important orientation(s) of linear structures in the image. Two test images along with their spectra are shown in Fig. 12. A spectrum σ λ (α) of an image f is computed as follows:

σ λ (α) = ∑ x∈Ω ( f ) γ α λ ( f ) (x). (6)
Again, the GPU implementation achieves a significant speedup (approximately 35 ×) for input images of 0.4 megapixels.

Rotation Detection of Music Sheet Scans

In the third practical application, we detected the rotation of music sheet scans. The music sheets were scanned and stored in an electronic archive. In the process of scanning, a random rotation could occur due to imperfect insertion of the paper to the scanner. The rotation was detected by computing linear closings with large SEs (λ = 250) of 81 different orientations within the angular range from -10 • to 10 • with the step of 0.25 • . The angular spectrum was computed according to Eq. ( 6). Two test images along with their spectra are shown in Fig. 13. In this case, the achieved speed-up was approximately 37 ×. It is known, that efficient GPU accelerations can only be obtained with simple, regular threads, using as low memory as possible. Hence, the separability principle is useful. Following this idea, one can form 2-D openings by concatenating 2-D erosion and 2-D dilation which are separable. A 1-D dilation algorithm with alike properties as the Bartovsky algorithm has been published in [START_REF] Dokládal | Computationally efficient, onepass algorithm for morphological filters[END_REF].

Assume a 2-D rectangular B in Eq. ( 1). It can be decomposed into two (horizontal and vertical) 1-D dilations and two 1-D erosions. Similarly, for hexagons we need to compute three erosions and three dilations, and for octagons four erosions and four dilations.

All these orthogonal operators need to be computed sequentially. One cannot expect to obtain the same performances in 2-D as with 1-D openings, since the execution times of the orthogonal operators are added together.

Conclusions

We have reviewed and compared the most efficient linear morphological opening/closing algorithms. At present, the fastest approaches (Van Droogenbroeck and Buckley, Morard et al., and Bartovsky et al.) compute the opening within a single image scan. The algorithm of Van Droogenbroeck and Buckley is the fastest one on CPU, however, it is efficient for 8-bit gray-scale images and for vertical and horizontal linear openings only. Morard and Bartovsky algorithms are applicable to any data accuracy (including floating point).

As described in the paper, both Morard and Bartovsky algorithms themselves are sequential. Hence, the only possibility of introducing more parallelism is on the thread execution level. We explain the choice of the algorithm (Bartovsky) to implement with regard to the GPU architecture (little memory, synchronous execution of threads). We have used various optimization techniques to speed up the code. Mapping various types of data (input, output and FIFOs) to various memory spaces is a crucial aspect. The choice of memory spaces is described in detail in Section 5.2. We also observed performance limits on small images, hence we try to introduce more spatial parallelism by splitting small images. Finally, we give rules of optimal configuration regarding input image size and the features of the available GPU. We provide comparisons with the fastest implementations on CPU and on GPU. The current state-of-the-art standard, OpenCV_GPU, is suitable for GPU architecture but has time complexity dependent on the SE size. The proposed implementation obtained stable performance over all orientations and sizes of the structuring element. For a single opening of a large image, we have measured up to 50 × speedup compared with CPU. For small images, the gain is significant (up to 37 × speedup) if one computes a set of openings in multiple directions. For example, within 60 ms, the GPU implementation is capable of computing a single opening at arbitrary angle of a 60 Mpx image, or a set of openings in 180 directions of a 640×640 px image.

To conclude, this solution is suitable for applications with large, industrial images, running under severe timing constraints, such as production control in e.g. metallurgy or textile industry. A typical such application requiring using high-resolution images, and running under severe time constraints is the surface control. Thin (often µm) cracks in large surfaces require using high resolution images, and the timing is given by the industrial cycle.

Source codes of the CPU and GPU implementations of the Bartovsky algorithm are publicly accessible under GNU GPL license [START_REF] Karas | CUDA-based Linear Openings[END_REF].

Fig. 1

 1 Fig. 1 Example application of linear openings. (a) input image, (b) linear opening with a digital line segment of length λ = 41 pixels and orientation 70 degrees, (c) enhancement of linear structures, and (d) color coded local orientation of linear structures calculated using definitions in Eq. (3). In (a) and (c), two zoomed regions are shown for a better comparison.

Fig. 2

 2 Fig. 2 Comparison of CPU implementations for horizontal and vertical SEs of size approx. 5% of the image width.

  [26], Soille et al. extended this work to arbitrary-oriented openings. In [27], Clienti et al. improved the HGW algorithm by removing the image backward scanning to reduce latency.

Fig. 3

 3 Fig. 3 Thread mapping and memory hierarchy in the GPU architecture.

Fig. 4 Fig. 5

 45 Fig. 4 The mapping of threads to the 2-D image grid in the GPU implementation of the algorithm for computing 1-D morphological openings/closings. Each thread, denoted by its ID, is mapped to an individual image column or row, respectively.

Fig. 6

 6 Fig. 6 Image split applied for opening/closing with vertical SE of size 2. Image is split into 2 zones, introducing twice more threads. The threads are denoted by vertical arrows, analogously to Fig. 4.

Fig. 7

 7 Fig. 7 FIFO storage in the shared and the global memory. The numbers indicate the order of pushing the element into the FIFO. The marked elements (the first one and the last two) are also cached in registers.

Fig. 8 Fig. 9

 89 Fig. 8 Comparison of optimized GPU implementations of Morard and Bartovsky algorithms. For comparison with the basic GPU implementation, refer to Fig. 5.

7

  Perspectives -Extention to 2-DThe 2-D opening is not separable into two orthogonal 1-D openings as is the dilation. Hence, one cannot directly combine two orthogonal 1-D openings to obtain a 2-D opening.

Fig. 12

 12 Fig.12Texture images[START_REF] Randen | Brodatz Textures[END_REF] and their angular spectra.

Table 2

 2 Comparison of CPU and GPU implementations for computation of an angular spectrum. The speedup is computed by comparing the GPU implementation with the best CPU (Morard) implementation.

	Image	Image size No. of		Time [ms]		Speedup
			angles CPU (Bart) CPU (Mor) GPU (Bart) GPU/CPU
	Fingerprint	784 × 1133	180	3768.7	3129.2	115.9	27.0
	D15	640 × 640	180	2425.8	2014.0	58.0	34.7
	D47	640 × 640	180	2397.7	2014.0	55.1	36.6
	Music1	1000 × 1411	81	5830.0	4840.7	129.3	37.4
	Music2	1000 × 1301	81	7297.1	4958.8	136.5	36.3

General-purpose computing on graphics processing units

available in Fulguro image processing library[START_REF] Ch | Fulguro image processing library[END_REF]

Acknowledgements

This work has been done during ERASMUS stage of Pavel Karas at ESIEE Paris. The presented work is the result of collaboration between ESIEE Paris (A3SI team), Mines-ParisTech (CMM) and Masaryk University in Brno (CBIA). This stay period of Pavel Karas at ESIEE Paris has been also supported by the Ministry of Education, Youth and Sport of the Czech Republic (Projects No. MSM-0021622419 and No. LC535). This research was also supported by the Grant Agency of the Czech Republic (Grant No. P302/12/G157).


G S Janniche

C Mouvet

H.-J Albrechtsen

VERTICAL SMALL SCALE VARIATIONS OF SORPTION AND MINERALIZATION OF THREE HERBICIDES IN SUBSURFACE LIMESTONE AND SANDY AQUIFER

Keywords: Isoproturon, mecoprop (MCPP), acetochlor, degradation, pesticide, fate

Vertical variation in sorption and mineralization potential of mecoprop (MCPP), isoproturon and acetochlor was investigated at low concentrations (µg-range) at the cm-scale in unsaturated sub-surface limestone samples and saturated sandy aquifer samples from an agricultural catchment in Brévilles, France. From two intact core drills, four heterogenic limestone sections were collected from 4.50-26.40 m below surface (mbs) and divided into 12 sub-samples of 8-25 cm length, and one sandy aquifer section from 19.20-19.53 m depth divided into 7 sub-samples of 4-5 cm length. In the sandy aquifer section acetochlor and isoproturon sorption increased substantially with depth; in average 78 % (acetochlor) and 61 % (isoproturon) per 5 cm. Also the number of acetochlor and isoproturon degraders (mostprobable-number) was higher in the bottom half of the aquifer section (93->16 000 /g) than in the upper half (4-71 /g). One 50 cm long limestone section with a distinct shift in color showed a clear shift in mineralization, number of degraders and sorption: In the two brown, uppermost samples, up to 31% mecoprop and up to 9% isoproturon was mineralized during 231 days, the numbers of mecoprop and isoproturon degraders were 1 300 to >16 000 /g, and the sorption of both isoproturon and acetochlor was more than three times higher, compared to the two deeper, grayish samples just below where mineralization (≤4 %) and numbers of degraders (1-520 /g) were low for all three herbicides. In both unsaturated limestone and sandy aquifer, variations and even distinct shifts in both mineralization, number of specific degraders and sorption were seen within just 4-15 cm of vertical distance. A simple conceptual model of herbicides leaching to groundwater through a 10 m unsaturated limestone was established, and calculations showed that a 30 cm active layer with the measured sorption and mineralization values hardly impacted the fate of the investigated herbicides, whereas a total thickness of layers of 1 m would substantially increase natural attenuation.

Introduction

Pesticides and their degradation products are often detected in surface water and groundwater [START_REF] Battaglin | Occurrence of sulfonylurea, sulfonamide, imidazolinone, and other herbicides in rivers, reservoirs and ground water in the Midwestern United States, 1998[END_REF][START_REF] Kolpin | Occurrence of selected herbicides and herbicide degradation products in Iowa's ground water[END_REF]Thorling, 2009) even above the EU groundwater guideline of 0.1µg/L (European Parliament and Council, 2006) and they are in this way posing a substantial risk to drinking water resources. This diffuse pesticide pollution has become one of the major environmental challenges for using groundwater as drinking water, and to predict the fate of pollutants in groundwater the entire catchment has to be considered [START_REF] Turner | Future trends in transport and fate of diffuse contaminants in catchments, with special emphasis on stable isotope applications[END_REF].

Knowledge on sorption and degradation is essential for evaluating risks of pesticides contamination of groundwater since these two processes are the major drivers in natural attenuation. Degradation of pesticides is often considered to decrease with depth [START_REF] Fomsgaard | Degradation of Pesticides in Subsurface Soils, Unsaturated Zone -A Review of Methods and Results[END_REF] but degradation can be faster in sub-soil than in top-soil [START_REF] Mills | Quantification of acetochlor degradation in the unsaturated zone using two novel in situ field techniques: comparisons with laboratory-generated data and implications for groundwater risk assessments[END_REF]. Sorption may also decrease with depth because of decreasing organic matter content which on the other hand may increase the bioavailability [START_REF] Bending | Microbial aspects of the interaction between soil depth and biodegradation of the herbicide isoproturon[END_REF]. Spatial variability in herbicide degradation and sorption has mainly been studied in topsoil or the upper meter of agricultural fields (e.g. [START_REF] Bending | Spatial heterogeneity in the metabolism and dynamics of isoproturon degrading microbial communities in soil[END_REF][START_REF] Bending | In-field spatial variability in the degradation of the phenyl-urea herbicide isoproturon is the result of interactions between degradative Sphingomonas spp. and soil pH[END_REF][START_REF] Charnay | Spatial variability in C-14-herbicide degradation in surface and subsurface soils[END_REF][START_REF] Jacques | Spatial variability of atrazine sorption parameters and other soil properties in a podzoluvisol[END_REF][START_REF] Fredslund | Spatial variation in 2-methyl-4-chlorophenoxyacetic acid mineralization and sorption in a sandy soil at field[END_REF][START_REF] Vinther | Field-scale variation in microbial activity and soil properties in relation to mineralization and sorption of pesticides in a sandy soil[END_REF][START_REF] Walker | Spatial variability in the degradation rate of isoproturon in soil[END_REF][START_REF] Walker | Spatial variability in the degradation rates of isoproturon and chlorotoluron in a clay soil[END_REF] and even at the millimeter scale [START_REF] Gonod | Spatial variability of 2,4-dichlorophenoxyacetic acid (2,4-D) mineralisation potential at a millimetre scale in soil[END_REF].

Although several studies have focused on pesticide fate in subsurface materials and aquifers only few have studied the vertical variability in these environments [START_REF] Albrechtsen | Degradation of herbicides in shallow Danish aquifers: an integrated laboratory and field study[END_REF][START_REF] Charnay | Spatial variability in C-14-herbicide degradation in surface and subsurface soils[END_REF][START_REF] Coquet | Pesticide adsorption in the vadose zone: a case study on Eocene and Quaternary materials in Northern France[END_REF][START_REF] Johnson | Potential for isoproturon, atrazine and mecoprop to be degraded within a chalk aquifer system[END_REF][START_REF] Larsen | Mecoprop, isoproturon, and atrazine in and above a sandy aquifer: Vertical distribution of mineralization potential[END_REF][START_REF] Wood | Spatial variability in herbicide degradation in the subsurface environment of a groundwater protection zone[END_REF]. Even fewer have considered chalk or limestone, probably because of the complexity and the difficulties in sampling such geological settings. Vertical variation in limestone of atrazine mineralization was investigated in seven samples collected from 4.6-10.0 mbs [START_REF] Charnay | Spatial variability in C-14-herbicide degradation in surface and subsurface soils[END_REF], and small-scale vertical variation in isoproturon degradation was investigated in chalk (four 5-cm core sections from 4.42-4-88 mbs) [START_REF] Johnson | Potential for isoproturon, atrazine and mecoprop to be degraded within a chalk aquifer system[END_REF], but since degradation was negligible in both studies, their contribution to insight in the spatial variation was limited. In limestone (Charnay) the sorption of atrazine, isoproturon and metamitron varied vertically only to a limited degree in 10 out of 11 samples (3-12 mbs), but one sample had significantly higher K d -values which could not be explained by the content of clay or organic matter [START_REF] Coquet | Pesticide adsorption in the vadose zone: a case study on Eocene and Quaternary materials in Northern France[END_REF]. Hence, only vertical variation in sorption but not in degradation of herbicides in subsurface carbonate rock was reported.

If sediment layers with increased mineralization could be identified it would be a step towards better understanding the interaction between sediment and biodegradation, which may improve the predictions and modeling of herbicides fate. Such sediment layers may in fact be "ecotones" where the microbial activity is high, found at the interface between different lithological layers [START_REF] Goldscheider | Review: Microbial biocenoses in pristine aquifers and an assessment of investigative methods[END_REF]. Spatial variability at the cm-scale in natural mineralization potential and sorption capacity of herbicides has to our knowledge not been studied before in deep subsurface samples (below 5 mbs), and especially not in limestone. The aim of this current work was to investigate i) if there is a vertical small scale spatial variability in sorption and in mineralization potential of the herbicides mecoprop, isoproturon and acetochlor in sediments from subsurface limestone and sandy aquifer; ii) if the variations could be related to a simple geochemical parameter; and iii) at which such a variation occurs and becomes important for the whole system.

Materials and Methods

Sediment sampling and characterization

Sediment was collected as intact cores in 2005 (Janniche et al., submitted), from a 4.5 km 2 agricultural catchment with diffuse pesticide pollution. The unsaturated zone consists of lutetian heterogeneous limestone, and the sandy aquifer is unconfined and aerobic with the groundwater table between 13.3 and 41.9 m below surface (mbs) [START_REF] Morvan | Pesticides in the groundwater of a spring draining a sandy aquifer: Temporal variability of concentrations and fluxes[END_REF]. Five sections were collected from two cores (Pz14 and Pz17c); four from the unsaturated limestone (sections A to D) and one from the saturated sand (section E). The sections were sampled from 4.50 to 26.40 mbs (Table 2). The sections were further divided vertically in 4 to 25 cm long subsamples. The four limestone sections were selected based on a visual inspection: three with rather consistent color and texture over depth, and one section with a distinct shift in color that divided the section in a brownish upper part and a grayish lower part. The sand section had a color gradient over the section length; ranging from brown in the upper samples to dark grey especially in the two bottom samples.

Water content of the sediment samples was measured as weight loss after 24 h at 105 °C. Sediment pH w was measured by electrode in 10 g sediment mixed with 10 mL deionized (Bigham et al., 1996). Total inorganic carbon (TIC) was calculated as the difference between total carbon content (TC) and Total organic (TOC). TOC and TC were measured based on the method by [START_REF] Heron | Determination of nonvolatile organic carbon in aquifer solids after carbonate removal by sulfurous acid[END_REF] by a total elemental carbon analyzer (LECO CS-225) in triplicates (0.2000 g) where carbonates were removed by adding 2 mL 6 % sulphurous acid (H 2 SO 3 ) six times to ensure that no bubble formation as released CO 2 could be observed.

Grain size (only for unconsolidated samples) was characterized by sieving (sieves of 0.063-2.0 mm) and, for the particles below 63 μm, by particle distribution (Sedigraph 5100, Micromeritics Gemini). The specific surface area was measured by Multipoint N 2 -BET analysis (Micromeritics, Gemini III 2375 surface area analyzer), after being outgassed (Micromeritics, FlowPrep 060 Degasser) for 4h at 70 °C.

Herbicides

Isoproturon and acetochlor are hydrophobic, non-polar organic herbicides whereas mecoprop is a hydrophilic phenoxy acid with a pK a value of 3.78 and at near-neutral pH it is therefore negatively charged (Table 1). The three herbicides have all been used in the field site to a wide extent. The experiments were carried out with [Ring-14 C]-pesticides (Table 1).

In the MPN also unlabelled acetochlor with 99.4 % purity was used.

Table 1

Sorption investigations

Sorption was investigated by a batch equilibrium technique based on the OECD guideline 106 (1993) at initial pesticide concentrations of 1 and 50 µg/L. Stock solutions of 14 C-labeled pesticide was prepared with 14 C-labelled pesticides in sterile filtered MilliQ water (Clausen et al., 2001). The sediment samples were crushed (only limestone), freeze-dried and sieved (<2 mm) [START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF]. The investigations were carried out in glass tubes (10 mL) with Teflon caps at 10°C in the dark and previous tests excluded adsorption on those materials [START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF]. The water from the Brévilles Spring in the catchment was used as equilibrium liquid in the investigations as representative for the groundwater after treatment with activated carbon to remove pesticides and sterile filtering. Triplicate 5 ± 0.0005 g of sediment was pre-equilibrated with 4 mL sterile filtered (0.2 µm Minisart, Sartorius) water for 20 h, then 1 mL pesticide solution was added resulting in an initial concentration of 1 or 50 µg/L (these low concentrations were chosen to mimic environmentally relevant concentrations), and each test tube was placed in a vertical rotator.

Based on preliminary investigations [START_REF] Mouvet | PEGASE. Pesticides in European Groundwaters: detailed study of representative quifers and Simulation of possible Evolution scenarios[END_REF] 168 h (for mecoprop and acetochlor) or 336 h (for isoproturon) was sufficient to reach equilibrium, and after this equilibrations time the suspension was centrifuged and 1 mL of supernatant collected for quantification of 14 C-activity (Clausen et al., 2001). Reference samples were prepared in triplicate without sediment but otherwise handled identically to the other samples. The amount of adsorbed pesticide was calculated as the difference between the 14 C-activity in the supernatant solution at equilibrium and the 14 C-activity in the reference solution (detection limit of 1 %). To check if the herbicides were degraded during the sorption equilibration 1 mL of supernatant was transferred to a "double-vial" system; 20 mL polyethylene scintillation vial containing a 6 mL scintillation vial with 1 mL 2.5 M NaOH to trap 14 CO 2 after acidification of the subsample (0.1 mL 37%). The inner-vial was removed after 48 hours, added scintillation cocktail and the 14 C-activity in e.g. the 14 CO 2 -fraction was quantified.

Mineralization

Mineralization was investigated in laboratory batch incubations set up under unsaturated or saturated conditions depending on the origin of each sample (Janniche et al., submitted) and incubated aerobically in the dark at 10°C for up to 251 days. Autoclaved (20 min at 1.5

bars and 120°C three times with one-day intervals) control incubations and duplicates were set up for selected sediments.

Unsaturated sediment was set up with 40 g (wet weight) sediment (limestone samples were "gently" crushed prior to a size that could pass the bottleneck with a diameter of 3 cm) in a 100 mL DURAN ® bottle with Teflon inlayer caps and spiked with 0.5 mL of pesticide stock solution to reach initial concentrations of 1 μg/kg. The produced 14 CO 2 was trapped in a test tube placed inside each bottle with 2 mL 0.5 M NaOH. When sampled, the sodium hydroxide was replaced and mixed with 10 mL scintillation fluid and the 14 C-activity quantified. The bottles were left open for 5 minutes in a laminar flow bench to replenish oxygen. Mineralization was calculated as the cumulated 14 CO 2 evolution over time as a percentage of initial 14 C-pesticide e.g. 14 CO 2 / 14 C 0 (%).

Aquifer sediment was set-up with 40 g (wet weight) sediment and 60 mL groundwater in a 118 mL glass serum bottle, sealed with a 1 cm thick butyl rubber stopper and crimp caps.

1 mL stock solution of 14 C labeled pesticide was added through the stopper resulting in an initial concentration in the water phase of approximately 1 μg/L of each pesticide. As for the sorption investigations, the groundwater was collected from the Brévilles spring, treated with activated carbon and sterile filtered. When sampling, 2 mL filtered (0.2 μm PTFE-filter;

Advantec/MFS 13HP) subsamples were acidified to strip off (from the water phase) 14 CO 2 that then was captured by NaOH in a"double-vial" system as described earlier. Pesticide content over time was calculated as the fraction of initial 14 C-pesticide corrected for pesticide removed during sampling. Mineralization over time was calculated as 14 CO 2 evolution (corrected for loss to headspace) as a percentage of the corrected initial 14 C-pesticide e.g.

14

CO 2 / 14 C 0 (%).

Enumeration of specific pesticide degraders

The specific pesticide degraders was enumerated by a modified (Janniche et al., submitted) most-probable-number method ( 14 C-MPN) [START_REF] Lehmicke | C-14-Most-Probable-Number Method for Enumeration of Active Heterotrophic Microorganisms in Natural-Waters[END_REF]. For each pesticide, sediment samples were mixed with pyrophosphate-buffer (1.2 mM Na 4 P 2 O 7 , pH = 7.0) and 10-fold serial dilutions (10 0 -10 -4 ) were prepared from which five replicates where made, resulting in 25 MPN-tubes to achieve 1 enumeration. The dilutions were added a complex mineral medium (ISO 7827, 1984) with pH = 6.8 and 25 µg/L 14 C-mecoprop or 14 Cisoproturon, or 10 µg/L 14 C-acetochlor and 20 µg/L acetochlor (due to late receipt of the acetochlor isotope the MPN-assay was set up with unlabeled acetochlor, and 3 month later was isotope added). The MPN-tubes were incubated 11-23 months in the dark at 15 ˚C and the 14 CO 2 evolved during incubation was trapped and measured by liquid scintillation counting. Controls were set up with sterile pyrophosphate-buffer without sediment suspension.

The MPN values and 95 % confides intervals were calculated by the computer software MPN calculator [START_REF] Curiale | MPN Calculater[END_REF] based on the number of positive scored tubes in each of the five dilutions. A positive score means that the MPN-tube contained bacteria capable of mineralizing the specific herbicide. The MPN-tubes were scored positive if the trapped 14 CO 2 was greater than the average of 10 sterile controls plus three times the standard deviation [START_REF] Lehmicke | C-14-Most-Probable-Number Method for Enumeration of Active Heterotrophic Microorganisms in Natural-Waters[END_REF]). The MPNs containing mecoprop were scored as positive if the radioactivity in the basetrap was >2.7% of the initially added amount of labeled; the values for isoproturon was 1.2% and 2.0% for acetochlor.

Results and Discussion

Sediment characterization

Core sections A-D were taken from the heterogeneous unsaturated limestone and section E from the sandy aquifer (Table 2). The heterogeneous limestone samples had relatively low values of total inorganic carbon (TIC): 17-69 % CaCO 3 with the lowest values found in the samples with marls. Though high TIC content may interfere with the total organic carbon (TOC) measuring, the results have high reproducibility, and are equivalent to reports for similar settings, e.g. 0.1-2.9 % TOC in limestone and 0-5.9 % in marlstone [START_REF] Coquet | Pesticide adsorption in the vadose zone: a case study on Eocene and Quaternary materials in Northern France[END_REF], who also reported CaCO 3 content of 69-88 % in limestone and 69-74 % in marlstone. The heterogeneous nature of our limestone may explain the relatively low TIC content compared to samples of pure limestone.

The limestone sections from Pz14 (A and C) had in general higher contents of water, TOC and TIC but a lower specific surface area than those from Pz17c (B and D). Besides differences observed between sections, differences were also observed at cm scale within a given section. Thus, limestone section D had a clear visual color shift and a clear shift in specific surface area with much higher specific surface area in the two uppermost samples (D1 and D2) compared to the two bottom samples (D3 and D4). Also the TIC content shifted with a lower TIC content in D1 and D2 compared to D3 and D4 expressing more pure limestone in D3-D4 than in D1-D2. Besides, TOC was much lower in the first sample, D1, than the other samples D2-D4 within this section.

Aquifer section (E) from Pz17c had high specific surface area but low TOC (Table 2).

The color gradually changed across the whole section from brownish to grayish and TIC slightly increased with depth; samples E1 and E2 had approximately half the TIC content of E3-E7.

In summary, the five core sections represented two major geological settings and different limestone settings; the limestone cores varied in textures, but the geology of the individual sections was very similar at the macro geological scale. Section D and E had clear small-scale variations at the micro scale in the sediment characterization parameters.

Table 2

Sorption

The sorption investigations included two concentrations (1 and 50 µg/L) and linear distribution coefficients (K d ) were calculated. The different concentrations only to a limited degree affected the K d -values (Fig. 1), and only noteworthy for isoproturon in the aquifer samples E3, E4 and E5, where sorption at the high concentration surprisingly was clearly higher (0.78-1.69 L/kg) than at the low concentration (0.52-0.87 L/kg). However, all isoproturon K d values from section E followed the same trend of increased sorption over depth (Fig. 1).

Fig. 1

In limestone samples, sorption was in general low for all three herbicides and so was the standard deviation reflecting high reproducibility (Fig. 1 A-D). In all limestone sections, isoproturon sorption was at or below the detection limit of 0.010 L/kg except in section D (K d : 0.03-0.25 L/kg). These findings are within the range reported for isoproturon sorption to limestone and chalk: <0.01-0.55 L/kg [START_REF] Besien | The transport and behaviour of isoproturon in unsaturated chalk cores[END_REF][START_REF] Coquet | Pesticide adsorption in the vadose zone: a case study on Eocene and Quaternary materials in Northern France[END_REF][START_REF] Johnson | Potential for aerobic isoproturon biodegradation and sorption in the unsaturated and saturated zones of a chalk aquifer[END_REF][START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF][START_REF] Mouvet | PEGASE. Pesticides in European Groundwaters: detailed study of representative quifers and Simulation of possible Evolution scenarios[END_REF]. The K d values for acetochlor in the limestone was <0.01-0.34 L/kg; literature K d values for acetochlor, in the range of 0.48-7.11 L/kg [START_REF] Hiller | Effect of soil and sediment composition on acetochlor sorption and desorption[END_REF][START_REF] Lengyel | Acetochlor as a soil pollutant[END_REF][START_REF] Ma | Field dissipation of acetochlor in two New Zealand soils at two application rates[END_REF], are predominantly reported for topsoil which explains their substantially higher K d values. The mecoprop sorption to limestone was also low (≤ 0.15 L/kg) and in the range (0-0.17 L/kg) reported in the literature for chalk and limestone (Clausen et al., 2001;[START_REF] Kristensen | Mineralization of aged atrazine and mecoprop in soil and aquifer chalk[END_REF][START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF].

The sorption of isoproturon and acetochlor to limestone was highest in D1 and D2 with higher specific surface area, and lower TOC than in the other limestone samples (Table 2). This indicated that isoproturon and acetochlor sorption in heterogenic limestone to a higher degree was controlled by mineralogy than by organic carbon content. Isoproturon sorbs to kaolinite, but not to quartz, calcite, aluminia (Clausen et al., 2001) or iron oxides (Clausen and Fabricius, 2001), and then sorption is strongly related to the content of organic matter and clay [START_REF] Coquet | Pesticide adsorption in the vadose zone: a case study on Eocene and Quaternary materials in Northern France[END_REF][START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF]. This supports that the increased isoproturon sorption was caused by the increased clay content within the limestone marls.

Sorption of acetochlor is likely to be controlled by the same parameters, but comparable literature data are limited for acetochlor sorption in subsurface sediments, low in TOC.

Sorption of mecoprop was not increased in the same sections as isoproturon and acetochlor because the sorption mechanisms of the anionic mecoprop differ from the nonionic herbicides isoproturon and acetochlor; mecoprop sorbs to positive sites on the mineral surface (Clausen et al., 2001).

In sandy aquifer samples sorption was in general substantial for acetochlor (0.25-5.12 L/kg) and isoproturon (<0.01-1.69 L/kg) except in E1 (Fig. 1), whereas mecoprop sorption was very low (≤0.12 L/kg) in the two investigated samples (E2 and E3). The sorption could not be explained solely by variations in the investigated geochemical parameters (clay content, organic and inorganic carbon content). The grayish color observed in the bottom two subsamples (E6 and E7) of the aquifer section indicated reduced conditions which may, however, explain the increased sorption capacity over depth because reduced organic matter induce a more hydrophobic surface whereto nonionic compounds better sorbs [START_REF] Clausen | Sorption of the herbicide dichlobenil and the metabolite 2,6-dichlorobenzamide on soils and aquifer sediments[END_REF]. K d values of 0-0.4 L/kg have been reported for mecoprop sorption onto sand [START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF][START_REF] Pedersen | Pesticide degradability in groundwater: Importance of redox conditions[END_REF][START_REF] Tuxen | Fate of seven pesticides in an aerobic aquifer studied in column experiments[END_REF], and our finding is within this range. K d values for isoproturon are reported to be 0.007-1 L/kg in sand [START_REF] Coquet | Pesticide adsorption in the vadose zone: a case study on Eocene and Quaternary materials in Northern France[END_REF][START_REF] Madsen | Pesticide sorption by low organic carbon sediments: a screening for seven herbicides[END_REF][START_REF] Pedersen | Pesticide degradability in groundwater: Importance of redox conditions[END_REF][START_REF] Tuxen | Fate of seven pesticides in an aerobic aquifer studied in column experiments[END_REF]; our results are in close agreement with these values.

Mineralization potential

Mecoprop, isoproturon and acetochlor degraders were detected in many samples (Fig. 2) which indicated a degradation potential of these compounds in the catchment. Furthermore, this suggested that the catchment was polluted with all three herbicides since another study only detected phenoxy acid (mecoprop, dichlorprop and 2,4-D) degraders in the part of the aquifer exposed to phenoxy acids (<40 µg/L) but not in the upstream pristine aquifer (de [START_REF] De Lipthay | In situ exposure to low herbicide concentrations affects microbial population composition and catabolic gene frequency in an aerobic shallow aquifer[END_REF]. The number of mecoprop degraders found in our investigation was at same level as observed by de Lipthay et al. ( 2003) but the number of isoproturon degraders was higher than reported by [START_REF] Bending | Microbial aspects of the interaction between soil depth and biodegradation of the herbicide isoproturon[END_REF].

Fig. 2

In all the 12 limestone samples mineralization of acetochlor was below 2%. Isoproturon was mineralized (up to 9%) in D1 and D2; in the remaining samples <1% was mineralized.

19-31 % mecoprop was mineralized in D1 and D2 and 4 % in D3 and C4 (average of two replicates with 6 and 2 %). The other samples exhibited mecoprop mineralization below 1 %.

The number of degraders was in general high in section D, but also in sample C3 the number of isoproturon degraders was extremely high (>16 000 /g). Some samples had high number of degraders without subsequent high mineralization (e.g. D1 (acetochlor), C3 (isoproturon)).

In the sandy aquifer (E), mecoprop was also mineralized substantially; 10-22% was recovered as 14 CO 2 with only 49-790 /g mecoprop degraders. The mineralization of both isoproturon (1-5%) and acetochlor (3-6%) was markedly lower than for mecoprop, despite higher numbers of isoproturon (27-2 200 /g) and acetochlor degraders (4-16 000 /g).

Small scale variability

The 19 samples from the 5 sections fell into 3 different cases: 1) In unsaturated limestone sections A, B and C sorption, mineralization and MPN were in general limited for all three studied herbicides (one exception was mecoprop sorption which was ≤ 0.01 L/k in sample A1 but 10-fold higher (0.11-0.15 L/kg) in A2 -the sample just below ). 2) Unsaturated limestone section D clearly showed small scale variations with a clear shift in sorption, mineralization and MPN. 3) Sandy aquifer section E had evident small scale variation in sorption but without a clear shift; it was rather a gradient, but also with variation in especially MPN.

Unsaturated limestone from section D

Sorption of isoproturon and acetochlor, and mineralization of mecoprop and acetochlor as well as number of mecoprop and acetochlor degraders were increased in the upper part of section D (Fig. 1 and2). A shift in color divided the section in two: the two brownish uppermost samples (D1 and D2) with high sorption (isoproturon: 0.17-0.25 L/kg and acetochlor: 0.22-0.34 L/kg), mineralization (mecoprop: 19-31 % and isoproturon: 9 %), and number of herbicide degraders (mecoprop: 1300-1400 /g and isoproturon: 4300->16000 /g) and on the other hand the grayish, lower samples (D3 and D4) with lower sorption (isoproturon: 0.03-0.07 L/kg and acetochlor: <0.01-0.11 L/kg), mineralization (mecoprop: 0.5-4 % and isoproturon: 0.2-0.5 %), and number of herbicide degraders (mecoprop: 1-6 /g and isoproturon: 3-520 /g). Already after 14 days of incubation the mineralization in D1 and D2 was noticeable higher than in D3 and D4 (Fig. 3). Furthermore, the 14 CO 2 evolution in D4

was as low as in the autoclaved incubations with sediment from D1-D4. Hence, a clear small scale shift in sorption and mineralization potential within only 25 cm was evident in the limestone.

Sandy aquifer section E

In the sandy aquifer samples sorption was low in the two uppermost samples, and increased gradually with depth. Acetochlor sorption increased in average 78 % pr every 5 cm in the sandy aquifer section (E). Also the isoproturon sorption increased gradually with depth (in average 61 % pr 5 cm) in the sandy aquifer (E) and hence the small scale variance in acetochlor and isoproturon sorption within these 33 cm was substantial. Despite the acetochlor mineralization in the sandy aquifer was rather constant, the number of acetochlor degraders increased dramatically in the three bottom samples; from 4-9/g in E1-E4 to 120-16 000/g in E5-E7. Also the number of isoproturon degraders was higher in the bottom samples (93-2 200/g in E4-E7) than in the top samples (27-71/g in E1-E3), and mineralization was faster in E1-E4 (2-5 %) than in E5-E7 (<1 %). For mecoprop however, the middle part (E3-E5) had the highest number of degraders; 490-790/g compared to 49-93/g, and mineralization

was also substantially faster in E5 (22 %) compared to the rest of the aquifer (10-15 %).

Mineralization of acetochlor and isoproturon in the saturated sandy aquifer samples (Fig. 2) could very well be negatively affected by the increased sorption with depth, leaving less herbicide available for mineralization, and could explain the lack of faster mineralization, despite very high number of degraders in some samples. This inverse relationship between sorption strength and mineralization was reported for MCPA [START_REF] Jensen | Sorption-controlled degradation kinetics of MCPA in soil[END_REF]. This could explain why isoproturon was mineralized (1.7-4.7 %) in the uppermost samples (E2-E4) but not in the bottom samples (<1 % in E5-E7). However, acetochlor which sorbed stronger did not follow that trend.

Conclusively, sorption and mineralization potential varied substantially within just a few cm in both limestone and sand. These small scale variations were found in sections with shifts in color, but could not be strongly related to the investiaged geochemical parameters.

Importance of an active limestone layer on groundwater contamination

The investigated four limestone sections revealed an overall picture of very low herbicide sorption and mineralization, except in the upper two 15 cm samples of section D,

where both sorption and mineralization potential was increased. This section contained marls (Table 2) and upper Lutetian marly layers are ubiquitous in the catchment, and may reduce the flow rate and subsequent solute transport [START_REF] Gutierrez | Long-term transfer of diffuse pollution at catchment scale: Respective roles of soil, and the unsaturated and saturated zones (Brevilles, France)[END_REF]. This marly layer with increased water content may be an environment favoring microbial activity.

To evaluate the relative impact on mecoprop, isoproturon and acetochlor of such an active layer with respect to sorption and mineralization in the unsaturated limestone a simple conceptual model was set up. Although both matrix flow and fracture flow may occur in heterogeneous unsaturated limestone matrix flow has in several studies shown to be predominant in deep groundwater sites and at steady-state flow [START_REF] Haria | Mechanisms of groundwater recharge and pesticide penetration to a chalk aquifer in southern England[END_REF][START_REF] Mathias | The significance of flow in the matrix of the Chalk unsaturated zone[END_REF][START_REF] Wellings | Recharge of the Upper Chalk Aquifer at A Site in Hampshire, England 1. Water-Balance and Unsaturated Flow[END_REF]. Slow transfer velocities of 0.5 to 1 m/yr have earlier been determined for our catchment based on profiles of water contents, tritium activities, and nitrate concentrations [START_REF] Gutierrez | Long-term transfer of diffuse pollution at catchment scale: Respective roles of soil, and the unsaturated and saturated zones (Brevilles, France)[END_REF]. We therefore assumed that the dominant flow occurred in matrix in the calculations of the conceptual model consisting of a 10 m unsaturated heterogeneous limestone profile with an active layer of a certain thickness (x a ).

Sorption and mineralization was confined to this active layer and hence microbial activity is assumed in the matrix. Such an active laver with sorption will cause retardation (R) of the herbicides which was calculated assuming matrix flow as: ( )

1 + = = =
1 - = R v
x Delay w a Mineralization of herbicide will reduce the risk of groundwater contamination in accordance with the observed mineralization patterns (Fig. 3) and the removal of herbicides due to mineralization was calculated based on first order mineralization kinetics as:

% 100 exp 1 ×                 - = w a v R x k removal
where k is the first order mineralization constant. K d -values and k-values are mean values as observed for D1 and D2 (Fig. 1 and2). Unfortunately, density (ρ) and porosity (ε) were not measured specifically for D1 and D2 due to lack of sample material but in 7 other samples from the catchment; ρ ranged from 2.69 to 2.73 kg/L and ε from 0.08 to 0.37 (Janniche et al., submitted). We have used the mean value of 0.18 for ε and 2.72 kg/L for ρ in this simple conceptual model calculation which results in a ρ b of 2.2 kg/L. If a lower porosity (ε) and hence a higher bulk density (ρ b ) was applied the retardation would have increased, and vice versa.

Three scenarios were set up: 1) a reference scenario with no active layer, 2) an active layer of 30 cm thickness, and 3) an active layer of 1 m thickness (Fig. 4). In the reference scenario the herbicides are transported at the same velocity as water (v w ), and with a water velocity of 1 m/yr it will take 10 years before breakthrough of the herbicides through the 10 m long profile. Presence of a 30 cm thick layer (scenario 2) will only to a very limited degree cause delay. It would take only 10.2 year for mecoprop to reach the groundwater (bottom of profile), 10.7 year for isoproturon and 11.0 year for acetochlor (Fig. 4). However, the mineralization in the layer would remove approximately 24 % mecoprop, 10 % isoproturon and 3 % acetochlor. The heterogeneous limestone could contain several active layers or one thicker layer. This may be realistic since the limestone settings in the field site are more than 30 m thick in 61 % and less than 10 m thick in 7 % of the surface of the catchment [START_REF] Morvan | Pesticides in the groundwater of a spring draining a sandy aquifer: Temporal variability of concentrations and fluxes[END_REF]. Assuming the summed extent of active layers to be 1 m over the whole 10 m profile (scenario 3), then the breakthrough times would be 10.6 years (mecoprop), 12.4 years (isoproturon) and 13.3 years (acetochlor) respectively, and the proportion of mineralized herbicide would then become 59 % mecoprop, 31 % isoproturon and only 9 % acetochlor (Fig. 4). In conclusion, if only one active layer of 30 cm thickness is present pr 10 m the retardation due to sorption will be negligible for all three herbicides and the removal will negligible for acetochlor but not for isoproturon and mecoprop, whereas the effect of both retardation and removal will be notable for all three herbicides if the layer is 3 times thicker.

Considering the thickness of unsaturated limestone in the investigated field site is up to 42 m thick it is not unlikely that the summed extend of an active layer will be greater than the observed 30 cm. Hence the importance of such active layers in the unsaturated limestone depends on their magnitude and consequently a very detailed knowledge of the limestone characteristics of the site is paramount.

Conclusions

Within few centimeters small scale vertical variability was present beyond measuring and sampling uncertainties showing clear shifts in sorption and microbial activity in both heterogeneous unsaturated limestone and sandy aquifer. None of the investigated parameters (TOC; TIC; specific surface area; sand, silt and clay content; water content; pH w ) individually could explain the variations. These variations followed a color shift or a color gradient in the core section.

The significance of these layers, or ecotones, with increased sorption and microbial mineralization depends on the specific herbicide and the extent of these layers, and it is therefore important to determine this. Such a detailed knowledge will improve the predictions of herbicide fate in the catchment if the geology is well described and included in reactive geochemical models.

Fig. 3

 3 Fig.3

  is the solute velocity [m/y], ρ b is the bulk density [kg/L] calculated as ρ b = (1ε)ρ, ε is the porosity [-], and K d is the linear sorption coefficient [L/kg]. The delay in break through due to sorption of the herbicides in the active layer was calculated as:

Fig. 4

 4 Fig.4

Acknowledgement

The European Union FP6 Integrated Project AquaTerra (Project no. GOCE 505428) under the thematic priority: sustainable development, global change and ecosystems funded this study together with the Technical University of Denmark (DTU). 14 C-labelled acetochlor was kindly provided by Dow Agroscience. The authors thank Ellinor Lindberg and Mona Refstrup for laboratory assistance, and Torben Dolin for the graphic assistance. Also thank to Bolette Nygaard and BRGM who took active part during drilling at Brévilles, and Alexis Gutierrez and Nicole Baran for always being helpful.