Judicaël Menant
email: jmenant@insa-rennes.fr

Muriel Pressigout

Luce Morin
email: lmorin@insa-rennes.fr

Jean-Francois Nezan
email: jnezan@insa-rennes.fr

Optimized fixed point implementation of a local stereo matching algorithm onto C66x DSP

Stereo matching techniques aim at reconstructing disparity maps from a pair of images. The use of stereo matching techniques in embedded systems is very challenging due to the complexity of the state-of-the-art algorithms. An efficient local stereo matching algorithm has been chosen from the literature and implemented on a c6678 DSP. Arithmetic simplifications such as approximation by piecewise linear functions and fixed point conversions are proposed. Thanks to factorisation and precomputing, the memory footprint is reduced by a factor 13 to fit on the memory footprint available on embedded systems. A 14.5 fps speed (factor 60 speed-up) has been reached with a small quality loss on the final disparity map.

I. INTRODUCTION

Embedded vision is the merging of two technologies corresponding to embedded systems and computer vision. An embedded system is any microprocessor-based system that is not a general-purpose computer [START_REF] Alliance | What is embedded vision[END_REF].

The goal of our work is to implement computer vision algorithms in modern embedded systems to provide them with stereo perception. Such systems take up little space and consume little power, so they are ideal for widespread integration into everyday objects. However, their architecture is significantly different to desktop systems, posing a challenge on how new algorithms and implementations that effectively and efficiently use these computational capabilities can be created. This paper focuses on Binocular Stereo-Vision algorithms. Stereo Matching aims to create 3D measurements from two 2D images, generating a disparity map which is inversely proportional to the distance of any object to the acquisition system. Such maps are used in scenarios where distance must be computed in all domains of computer vision. This is a strategically important knowledge field to test computation acceleration with.

Active devices such as Kinect [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF] are able to produce real time depth map. This kind of device works by emitting an infra-red grid on the observed scene. The depth map is then deduced from this sensed-back grid. Those devices are limited to indoor use with a 5 meter range [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF] and they are sensitive to infra-red interferences. Binocular stereo matching algorithms bypass these limitations and they are multi-purpose.

Energy-efficient embedded platforms are now available. The c6678 platform is an 8-core DSP platform clocked at 1 GHz with a standard 10W power consumption. The challenge is now to find and adapt algorithms and implementations that can fully exploit the powerful computational capabilities of such an architecture.

Most existing implementations of stereo matching algorithms are carried out on desktop GPUs, leading to poor energy efficiency. To be efficiently implemented on embedded systems, algorithms must be ported to fixed point implementation. Fixed point implementations lead to quantization noise and quality loss, thus a trade-off between precision and quality must be found.

In this paper, an efficient fixed point implementation of a stereo matching algorithm on the c66x architecture is proposed. First, the c66x core and its architecture are introduced in section II. Then the studied stereo matching algorithm is exposed in section III. The studied algorithm is a local search modified in order to fit the c66x architecture. Next the quality and speed results are presented in section IV. Finally the perspectives and future work conclude the paper in section V.

II. THE C66X MULTI-CORE DSP PLATFORM

The c6678 platform is composed of 8 c66x DSP cores and it is designed for image processing.

A. Memory architecture

Memory is a critical point parameter in embedded systems. Large memories being slower than smaller ones, modern systems integrate several memory layers in order to increase the memory capacity without increasing the access time to data. The memory hierarchy of the c6678 platform is exposed described below :

• 512 MBytes of external DDR3 memory. This memory is slow, and the bus bandwidth is limited to 10 GBytes/s. This memory is shared between all 8 cores. It is connected to the 64-bit DDR3 EMIF bus. In most digital signal processing applications, a loop kernel is a succession of interdependent sequential operations. These loops must be pipelined to be parallelized.

Loop pipelining is done by the compiler, and a factor 6 speed-up can be easily achieved with a little human work. However a pipelined loop can not have any jumps, and conditional statements (if/else) must be avoided inside the loop. These rules must be kept in mind when writing efficient code and designing algorithms.

2) SIMD: Single Instruction Multiple Data (SIMD) are instructions that are executed on multiple data. A SIMD instruction considers one or two registers (32 or 64 bits respectively) as a set of smaller words. For instance a 32bit register can be seen as a group of 4 8-bit words, and is called a 4-way 8-bit SIMD instruction.

This kind of instruction is very useful in image processing because most of the manipulated data are 8-bit (pixels).

C. FPU

The FPU (Floating Point Unit) is a logic unit which is able to execute operations on floating point numbers. The c66x has a basic FPU ; Nevertheless, this FPU has low support of SIMD (two ways maximum), whereas SIMD on fixed point numbers is up to 8 ways. Moreover a floating point number is always 4 byte wide, and causes a higher memory usage.

III. PROPOSED ALGORITHM

A stereo matching algorithm computes depth information from two cameras. The goal of these algorithms is to match a pixel in left image with one in the right image. Due to their high computational complexity, designing real time stereo matching algorithms is a great challenge.

In the past decades a lot of work has been done on stereo matching algorithms to increase their quality and time efficiency. There are two main categories of stereo matching algorithms : local and global methods [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF].

Global methods minimize global energy on the entire disparity map or scanline. Global methods are very good in terms of quality, but not very well suited for real time applications.

Local methods are area-based algorithms, they match pixels by taking into account the neighbourhood of these pixels. Local methods give less accurate results than global methods, but they are much faster.

In this paper, only local methods are considered for real time considerations. The Middlebury Stereo Vision Website [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] is a reference for the comparison of Stereo Matching algorithms. The best ranked algorithms of this database and the related papers has been studied. As the highest ranked algorithm is not published, the algorithm proposed by Mei [START_REF] Mei | On building an accurate stereo matching system on graphics hardware[END_REF] has been targeted for a DSP Implementation.

First, the original algorithm [START_REF] Mei | On building an accurate stereo matching system on graphics hardware[END_REF] will be introduced, and all algorithm blocks will be detailed. Then, real time limitations of this algorithm will be discussed and solutions to bypass those limits will be exposed. To finish this section, some implementation details will be given.

A. Algorithm overview

Figure 1 represents the chosen algorithm [START_REF] Mei | On building an accurate stereo matching system on graphics hardware[END_REF]. It is composed of three main steps :

• The cost construction which computes matching cost for each possible disparity and for all pixels. • The cost aggregation which refines those cost maps thanks to an iterative algorithm. • The disparity selection which selects the minimum matching cost (i.e. the best match) and deduces the disparity.

A median filter is applied at the output of the disparity selection step. The median filter is common and will not considered in this paper.

1) Cost construction:

The cost construction step takes left and right images and computes a matching cost for all possible disparity levels. Its output is a cost map for each disparity level (same size as input image). It is computed thanks to equation [START_REF] Alliance | What is embedded vision[END_REF].

For a disparity d, a pixel p of coordinates (x, y) is taken in the left image and compared to a pixel p d of coordinates (x + d, y) in the right image. The result of the comparison is an error that is called matching cost. Equation (1) describes the matching cost that is computed on each pixel and for all possible disparities.

Cost(p, d) = C T AD (p, d).λ T AD + C CEN (p, d).λ CEN (1)
Equation ( 1) is composed of two sub-costs summed together : Fig. 2: 8 bits census signature example a) Truncated absolute difference: The C T AD expresses the similarity of two pixels, it focuses only on grey levels similarity, and not about its neighbourhood.

C T AD (p d , d) = 1 -e -|I l (p)-Ir (p d )| thr (2)
I r and I l are left and right grey level images. thr is a threshold (i.e. a constant) which defines the useful dynamics of truncated absolute difference. It must be superior to the noise value in image, 20 is a common value in the literature [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF].

b) Census Cost: The census gives an information about local texture that C T AD does not provide.

Census produces an 8-bit signature for each pixel of an input image. As shown in figure 2, this signature is obtained by comparing each pixel to its 8 neighbours. The census is referred as cen l , cen r in equation ( 3) for respectively left and right grey level images.

C CEN (p, d) = 1 8 7 k=0 0 if cen r (p)[k] = cen l (p d )[k] 1 if cen r (p)[k] = cen l (p d )[k] (3)
In equation ( 3), C CEN is the cost associated to census. cen l (p)[k] and cen r (p)[k] are the k th bits of the 8-bit census signature for pixel p in the left and right images respectively. p d is the pixel of coordinates p -d, with d the disparity level.

Equation ( 3) is a sum of bit-to-bit comparisons, and these bit-to-bit comparisons are equivalent to a boolean operation called "exclusive or" (XOR).

The output of the cost construction step is one cost map per disparity level. Those cost maps are the input of the cost aggregation step.

2) Cost aggregation: The cost construction step has a low computing cost, but it provides noisy matching cost maps. This noise is mainly produced by C CEN which is random when compared regions are not correlated. To remove this noise, the cost aggregation step performs smoothing on areas with similar colour in the original image. The cost aggregation step is performed independently on each cost map. This is a key point regarding implementation.

The cost aggregation algorithm's structure is similar to a bilateral filter. The cost aggregation step is performed iteratively with varying parameters, it is defined by equation ( 4).

E i+1 (p) = W (p, p + )E i (p + ) + E i (p) + W (p, p -)E i (p -) W (p, p + ) + 1 + W (p, p -) (4)
E i is the current cost map to be refined, E 0 is the output of cost construction (equation ( 1)).

Pixels p + and p -have a position relative to pixel p : 4) is computed alternatively for horizontal and vertical aggregation :

• p + = p + ∆ i • p -= p -∆ i Equation (
• When i is odd, it is a vertical aggregation, the offset ∆ i is vertical. • When i is even, it is a horizontal aggregation, the offset ∆ i is horizontal.
At each iteration the parameter ∆ i grows, thus further pixels p + and p -are used for smoothing p. ∆ i evolves according to equation ( 5), the influence range is limited by the modulo (here ∆ i ∈ [0, 32]).

∆ i = f loor(i/2) 2 mod 33 (5)
Weights W in equation ( 4) are defined by equation ( 6) :

W (p 1 , p 2 ) = e C d .∆i-sim(p 1 ,p 2 ) L2 (6) sim(p 1 , p 2 ) = col∈(r,g,b) (I r col(p 1 ) -I r col(p 2 )) 2 (7)
In equation ( 6), C d is a weight applied to distance [START_REF] Mei | On building an accurate stereo matching system on graphics hardware[END_REF] and L2 is the weight applied to similarity [START_REF] Mei | On building an accurate stereo matching system on graphics hardware[END_REF]. I r {r, g, b} and I l {r, g, b} are the RGB (Red, Green, Blue) signals of right and left images.

3) Disparity selection: The disparity selection step minimizes the matching cost. To do so, the Winner Takes All (WTA) strategy [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] is applied. The WTA strategy is a simple arithmetic comparison expressed by equation (8).

Disp(p) = argmin d∈[0,N disp] E d,N it (p) (8)
The output of disparity selection is a dense integer disparity map providing a disparity value for each pixel in the right image.

B. Arithmetic simplification

Section III-A has exposed the original algorithm. Now the proposed modifications in order to speed up the algorithm without too much degradation are exposed.

To be fast, all algorithm blocks must be implemented in fixed-point arithmetic, because the DSP is more efficient with fixed points (see II-C).

To be easily implemented in fixed-point arithmetic all functions must be described with classical operators (+, -, * ). Functions such as square root or exponential must be avoided. These principles are applied in the proposed modifications below.

1) Cost construction: Section III-D2 will show that C CEN fits well on a DSP architecture. But C T AD should be simplified (exponential function is not recommended). A common way in the literature [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] is to use a raw saturation described by equation (9).

C T AD (i, j, d) = min(thr, | I ri (i, j) -I le (i -d, j) |) (9)
Equation ( 9) is used instead of equation [START_REF] Khoshelham | Accuracy and resolution of kinect depth data for indoor mapping applications[END_REF].

2) Cost aggregation: Cost aggregation takes up 85 % of execution time (weight and aggregation rows in table III). Optimization of cost aggregation is thus most important in the stereo matching algorithm.

a) Number of iterations: The cost aggregation execution time being linear with respect to the number of iterations, a way to reduce execution time is to reduce the number of iterations. This point will be explored in section II.

b) Weight computing: Weight computing is quite complex because of the exponential function and the square root in equation [START_REF] Mercat | Implementation of a stereo matching algorithm onto a manycore embedded system[END_REF]. In this equation, the weights increase when two pixels are similar and decrease with distance. Equation (10) proposes a piecewise linear function to approximate equation [START_REF] Mercat | Implementation of a stereo matching algorithm onto a manycore embedded system[END_REF].

W (p 1 , p 2 ) = thr -trunc thr (sim(p1, p2)) thr .(1-∆.Cd) (10)
Equation ( 11) avoids use of square root by using a norm L 1 instead of a norm L 2 in the original equation [START_REF] Zhang | Real-time gpu-based local stereo matching method[END_REF].

sim(p1, p2) = col∈(r,g,b) | I r col(p 1 ) -I r col(p 2 ) | (11)
In equation (10), the denominator of the division is a constant and thus can be replaced by a multiplication by 1/thr. c) Aggregation: The core of cost aggregation is a simple weighted sum (numerator of equation ( 4)). The denominator of equation ( 4) expresses a normalization.

The normalization is important because it prevents divergence and keeps a constant average cost value. A constant average value is necessary because of the lack of signal dynamic due to fixed point. To reach this goal, the sum of weights must be equal to one.

This normalization is a critical point because it implies a division which is not handled by the hardware.

The proposed approximation to normalization is described by equation ( 12). The wanted property, the fact that the sum of weights equals one, is verified by equation (12

) (W p + W m + W o = 1). W p = W (p, p + ) 4 W m = W (p, p -) 4 (12a) W o = 1 -(W p + W m) (12b) E it+1 (p) = W p.E it (p + ) + W o.E it (p) + W m.E it (p -) (12c)
Equation ( 12) uses only basic operators : +, -, * and ≫ (shift right) for the division by 4.

C. Memory optimization

The main difficulty of implementation in embedded systems is the efficient use of the memory architecture. Data in the slow external memory (DDR3) should be avoided (see II-A). Best performances are reached when the memory footprint is lower than the internal memory size in order to bypass the external memory. The memory footprint of the algorithm has to be reduced. Next sections describe the proposed strategies to reduce memory footprint.

1) Disparity selection: Cost construction and cost aggregation are done independently for each disparity level (equations ( 1) and ( 4)). Winner Takes All (WTA) strategy (equation ( 8)) finds the minimum value of their outputs.

It is not necessary to store each level of disparity in memory when implementing argmin on the fly. Only the cost of the current disparity level (E d ), the current argmin d and minimum values (min d ) have to be kept in memory.

2) Weight precomputing: The iterative equation ( 12) is the core of this algorithm. Weights are independent of disparity but they are dependent on ∆ i (see equation ( 10)). That is why weights can be precomputed and reused for each disparity level.

Moreover, memory and computational complexity can be avoided. Indeed equation 11 implies W (p1, p2) = W (p2, p1) which also implies that W (p, p + ) and W (p, p -) are identical with an offset of ∆ i . This can easily be proved by taking a second point p ′ in weights table that has an offset of ∆ i :

p ′ = p + ∆ i = p + p ′ -= p ′ -∆ i = p W (p, p + ) = W (p + , p) = W (p ′ , p ′ -)
Because the tables W (p, p + ) and W (p, p -) are identical with an offset of ∆ i , only one of them needs to be computed. This property decreases memory requirement and computation by a factor two when precomputing weights.

Two weights tables are required per value of ∆ i (one horizontal and one vertical). Using this precomputing, each pixel operation is reduced to three additions and two multiplications plus normalization (see III-A2).

Finally all memory optimizations lead in a 13 times smaller memory footprint (for a 20 disparity level image) :

• A factor 4 is obtained by replacing floating point by fixed point. • Precomputed weights take 2 buffers per iteration level instead of 6. • It is no more necessary to store each disparity level in a buffer.

D. Hardware optimization

Previous parts were dealing with generic optimizations ;in this part, c66x-specific hardware optimizations are introduced.

As explained in section III-B, a set of basic operators must be used (+, -, * ). Those operators are available on all platforms, and often in SIMD version. But the c66x DSP has special instructions that are very useful for this application.

The following sections describe how to use these instructions efficiently.

1) Truncation: In part III-B, some truncation (minimum) functions were used (in equations ( 9) and ( 10)). The DSP has a 4-way 8-bit SIMD minimum instruction (minu4 [START_REF]TMS320C66x DSP, CPU and Instruction Set, Texas Instruments[END_REF]), so any saturation function can be implemented very efficiently without any conditional statement, thus this implementation is compatible with loop pipeline (see II-B1).

2) Census: The census cost is implemented very efficiently using special instructions.

When looking at the description of the census cost in section III-A1b, it cannot be implemented with usual operators, and might look hard to implement. But there is a very efficient implementation on the c66x DSP.

First, comparisons with a pixel's neighbourhood are not a problem, because they are precomputed once, and the cumulated execution time is negligible.

Second, an exclusive or can be used to compare the two 8-bit signatures in Equation (3). To count the number of bits that are equal to one, the "bit count" instruction can be used efficiently. This instruction takes a word and returns the number of bit which are equal to one. This instruction is available on 4-way 8-bit SIMD (bitc4 [START_REF]TMS320C66x DSP, CPU and Instruction Set, Texas Instruments[END_REF]).

IV. RESULTS

This section exposes the results of the proposed modifications in quality and performance (execution time). The fixed point implementation introduces more computing noise, thus the quality loss due to fixed point must be quantified.

A. Quality Assessment

The quality of algorithm output has been evaluated with Middlebury evaluation tools [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. The Middlebury evaluation gives the number of bad pixels in different parts of the image. In this paper two values are considered for quality assessment:

• The global number of bad pixels. It gives a good representation of the overall image quality. • The number of bad pixels in non-occulted areas of the image. This algorithm does not deal with occulted areas, and the disparity map on those areas is very noisy, thus results on occulted areas are partially random. Comparing random results together does not give clear results. That is why, when comparing two versions of this algorithm, it is better to look at non-occulted parts of images. 1) Data Set: Seven images from the Middlebury database [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF] are used to evaluate quality for the proposed algorithm (Aloe, Art, Laundry, Map, Sawtooth, Tsukuba, Venus).

2) Algorithm parameters: The algorithm's default constants from the literature are used [START_REF] Mercat | Implementation of a stereo matching algorithm onto a manycore embedded system[END_REF], [START_REF] Zhang | Real-time gpu-based local stereo matching method[END_REF], [START_REF] Daniel Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. These constants have not been tuned to increase quality on this particular data-set :

• thr in equation ( 9

B. Impact of approximations on quality

Impact of approximations described in part III-B on quality are exposed in table I.

The impact of weights and cost construction approximation on quality is 0.04 % on non-occulted pixels, and it can be considered as negligible. This can be explained by the fact that the equation (11) describes correctly the similarity of two pixels. Moreover, because outputs of those functions are only used for arithmetic comparisons in disparity selection (see equation ( 8)), impact on quality is minor. Since functions are monotonic (if arithmetic precision is not taken into account).

The normalization approximation is very intrusive, the quality loss is expected to be important on this point. Indeed the normalization approximation (equation ( 12)) degrades the quality by 1 % (see table I) on not occulted pixels, this is an acceptable quality lost.

C. Cost aggregation iterations

As explained in part III-B2, a easy way to increase performance is to reduce the number of iterations. Table II exposes the impact of number of iterations and the value of ∆ i on quality. The more iterations there are the better the quality is. But with more than four iterations, the quality does not increase very fast. This can be explained by the iterative algorithm which makes the result tend to an optimal value.

The execution time is almost linear with respect to the number of iterations. The sixth line of table II is a good tradeoff between the number of iterations and quality.

D. Fixed point precision

When implementing on the DSP using SIMD, some significant bits are removed to increase speed. Table IV exposes the quality lost when porting on DSP.

Results of only one image on DSP is presented, the quality loss is the same for others images. Indeed the quality loss is only due to precision lost on the DSP, and the loss is the same magnitude for any image. Results of image sawtooth are exposed in table IV.

The quality degradation from line 2 to 3 of table IV is caused by the reduction of arithmetic precision to enhance the speed of SIMD instructions. When using 8-bit data instead-of 16-bit data, twice as many operations can be achieved in one SIMD cycle.

Table IV shows a 0.5 % quality degradation caused by the loss of arithmetic precision. This quality degradation is not prohibitive. The stereo pair and output disparity map are shown figure 3.

E. Speed

In previous sections, the quality lost due to optimizations and arithmetic simplifications has been exposed. This last section deals with the speed performances resulting from those simplifications.

Table III shows execution time of each algorithm step, and the percentage of time used by these steps. The original version of the algorithm (floating point) runs at 0.24 fps on DSP, the fully optimized version runs at 14.5 fps, the speed-up factor is 60. Disparity selection step uses inherently fixed point (deals with integer disparity value), thus it is tagged as not implemented in fixed-point column. Disparity selection step can not benefit of loop pipeline without SIMD implementation, thus it is tagged as not implemented in pipeline column. Some steps were not optimized with SIMD because their execution times are too small to have a significant impact on total execution time.

Proposed modifications in section III are introduced at fixed-point column of table III. However, speeds of pipelined and SIMD versions could not be reached without proposed optimizations. For instance, as explained in section III-B2c, the cost aggregation core could not be pipelined because of normalization. That is why the speed-up of proposed modifications is not reduced to the speed-up between the original algorithm and its fixed point version.

To conclude this section on performance, the proposed modifications lead to a factor 60 speed-up (table III) with a quality of 8.68 % (table I) of bad pixel on image compared to 7.01 % before modification.

V. CONCLUSION AND PERSPECTIVES

The goal of this paper was to prove the feasibility of realtime dense stereo matching on a DSP platform. Currently, a cost-aggregation based algorithm is running at 14.5 fps on one c66x core. This speed has been reached with minor quality loss.

If the quality degradation is too important for applications, a trade-off depending on application can be found to increase quality. More significant bits can be allocated for computing with a quality improvement, but less operations can be done simultaneously, and eventually, more memory would be used (if intermediate data are store to 16 bits instead of 8). Moreover an easy way to increase the quality is to increase the number of iterations in the cost aggregation core.

The current implementation of the algorithm described in this paper uses only one of the eight cores of the c6678 platform. Future work will consist in using all the cores of this platform.

Fig. 1 :

 1 Fig. 1: Description of the Stereo Matching Algorithm

  ) : 20 • thr in equation (10) : 60 • Cd in equation (10) : 0,015 • λ cen in equation (1): 0,55 • λ tad in equation (1) : 0,45

Fig. 3 :

 3 Fig. 3: Inputs and output of the Stereo Matching algorithm

TABLE I :

 I Quality of algorithms

	Algorithm	All	Bad pixels Not occulted
	original	7.01 %	4.68 %
	weight	6.70 %	4.72 %
	weight and normalization	8.68 %	5.69 %

TABLE II :

 II Influence of iterations on cost aggregation.

	Iteration steps (∆ i )	Iterations	All	Bad pixels Not occulted
	1, 4, 9, 25	4	7.6 %	5.1 %
	1, 4, 9, 16, 1	5	7.7 %	5.3 %
	1, 4, 9, 16, 25	5	7.1 %	4.9 %
	1, 4, 9, 16, 25, 3	6	6.7 %	4.7 %
	1, 4, 9, 16, 25, 36	6	6.6 %	4.7 %
	1, 2, 4, 8, 16	5	7.8 %	5.4 %
	1, 2, 4, 8, 16, 32	6	6.7 %	4.6 %

TABLE III :

 III Execution times on DSP.

	algorithm step	execution number	floating point	fixed-point	pipeline	SIMD
	aggregation	N it.N Disp (95)	13,24 30,13 %	8,83 75,80 % 1,44 63,58 % 0,41 56,72 %
	cost construction	N Disp (19)	100,09	45,56 %	5,44 9,34 %	1,89 16,72 % 0,63 17,40 %
	disparity	N Disp -1 (18)	1,56 0,23 %	-	2,54 %	-	13,04 % 0,20 5,12 %
	weight	N it (5)	199,73	23,92 % 19,37 8,75 %	1,93 4,48 %	-	13,97 %
	total (FPS)			0,24		0,9		4,7 *	14,5
	* : data not directly measured (computed from measured values). -: not implemented parts. Units are machine million cycles
	(The DSP is clocked 1 GHz, 1 million cycles = 1 millisecond).

TABLE IV :

 IV Sawtooth quality

	Algorithm	All	Bad pixels Not occulted
	original	2.43 %	1.12 %
	normalization	3.41 %	1.35 %
	DSP (SIMD)	4.02 %	1.68 %

• The truncated absolute difference cost (C T AD ) expresses the similarity of two pixels and is described in section III-A1a. λ T AD is a weight associated to C T AD cost. • The census associated cost (C CEN ) gives an information about local texture and is described in section III-A1b. λ CEN is a weight associated to C CEN cost.


Majken C Haarder

Looms

D Jougnot
email: damien.jougnot@unil.ch

Niklas Linde
email: niklas.linde@unil.chelineb.

E B Haarder

M C Looms

Monitoring of saline tracer movement with vertically distributed self-potential measurements at the HOBE agricultural test site, Voulund, Denmark

Monitoring of saline tracer movement with vertically distributed self-potential measurements at the HOBE agricultural test site

Introduction

Quantification of water fluxes in the vadose zone is essential for many hydrological and environmental applications. Classical approaches based on matric potentials or tracer test data (e.g., [START_REF] Vereecken | On the value of soil moisture measurements in vadose zone hydrology: A review[END_REF][START_REF] Tarantino | Field Measurement of Suction, Water Content, and Water Permeability[END_REF] are limited by the punctual nature of such measurements. The data might be strongly influenced by local heterogeneities and water fluxes are obtained indirectly by differencing. [START_REF] Lunati | On the use of spatially discrete data to compute energy and mass balance[END_REF] demonstrated that large errors occur when mass and energy balances are computed from discrete measurements. One way to overcome the influence of local heterogeneity is to use geophysical measurements that are representative of larger volumes (e.g., [START_REF] Hubbard | Treatise on Water Science[END_REF]. The present contribution focuses on the self-potential (SP) method and to what extent it can be used to infer water fluxes and tracer transport at an experimental research site. The self-potential method is non-invasive and sensitive to subsurface flow and transport processes (for reviews, see [START_REF] Jouniaux | Review of self-potential methods in hydrogeophysics[END_REF][START_REF] Revil | Review: Some low-frequency electrical methods for subsurface characterization and monitoring in hydrogeology[END_REF][START_REF] Revil | The Self-Potential Method: Theory and Applications in Environmental Geosciences[END_REF]. It is a passive method, in which spatial and temporal variations of the electrical potential field are measured with respect to a reference electrode. The recorded self-potential data are given as a superposition of several contributions, which makes interpretation challenging.

The electrokinetic (EK) contribution (often referred to as the streaming potential) is directly related to the water flux and the properties of the electrical double layer found at the mineralpore water interface. Water flowing through the pore drags a fraction of the excess charge, which gives rise to a streaming current and a resulting electrical potential field. Electrokinetic effects have been studied for more than a century [START_REF] Helmholtz | Studien über electrische Grenzschichten[END_REF] and are well understood in water saturated porous media (e.g., [START_REF] Jouniaux | Review of self-potential methods in hydrogeophysics[END_REF][START_REF] Revil | The Self-Potential Method: Theory and Applications in Environmental Geosciences[END_REF]. Two main approaches have been proposed to simulate streaming current generation at partial saturation. The first focuses on how the streaming potential coupling coefficient varies as a function of water saturation (e.g., [START_REF] Guichet | Streaming potential of a sand column in partial saturation conditions[END_REF][START_REF] Revil | Streaming potentials in two-phase flow conditions[END_REF][START_REF] Darnet | Modelling streaming potential (SP) signals induced by water movement in the vadose zone[END_REF][START_REF] Allègre | Streaming potential dependence on water-content in Fontainebleau sand[END_REF], while the second focuses on how the excess charge dragged by the water varies with water saturation (e.g., Linde et al., 2007a;[START_REF] Revil | Electrokinetic coupling in unsaturated porous media[END_REF][START_REF] Linde | Comment on "Characterization of multiphase electrokinetic coupling using a bundle of capillary tubes model[END_REF][START_REF] Jackson | Multiphase electrokinetic coupling: Insights into the impact of fluid and charge distribution at the pore scale from a bundle of capillary tubes model[END_REF][START_REF] Mboh | Coupled hydrogeophysical inversion of streaming potential signals for unsaturated soil hydraulic properties[END_REF][START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF]. The lack of agreement between researchers is partly due to a limited number of well-controlled experiments and the multiple contributions to the measured signal, including electrode effects (e.g., [START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF].

Most well-controlled laboratory studies have focused on drainage experiments (Linde et al., 2007a;[START_REF] Allègre | Streaming potential dependence on water-content in Fontainebleau sand[END_REF][START_REF] Mboh | Coupled hydrogeophysical inversion of streaming potential signals for unsaturated soil hydraulic properties[END_REF] or drainage-imbibition cycles [START_REF] Haarder | Accepted. Estimation of recharge from long-term monitoring of saline tracer transport using electrical resistivity tomography[END_REF][START_REF] Vinogradov | Multiphase streaming potential in sandstones saturated with gas/brine and oil/brine during drainage and imbibition[END_REF][START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF][START_REF] Allègre | Self-potential response to drainage-imbibition cycles[END_REF]. Due to experimental difficulties, only few data sets describe the streaming potential coupling coefficient at partially saturated conditions at steady state [START_REF] Guichet | Streaming potential of a sand column in partial saturation conditions[END_REF][START_REF] Revil | Streaming potentials in two-phase flow conditions[END_REF][START_REF] Revil | Electrokinetic coupling in unsaturated porous media[END_REF]. These measurements suggest that the dependency on the streaming potential coupling coefficient with saturation is media dependent. [START_REF] Doussan | Variations of self-potential and unsaturated water flow with time in sandy loam and clay loam soils[END_REF] instrumented a lysimeter and monitored SP data and vertical water flux in partially saturated conditions. A limited number of hydrological studies have used SP monitoring data conducted at the field scale (e.g., [START_REF] Thony | Field characterization of the relationship between electrical potential gradients and soil water flux[END_REF][START_REF] Perrier | Characterization of electrical daily variations induced by capillary flow in the non-saturated zone[END_REF][START_REF] Revil | Self-potential signals associated with variations of the hydraulic head during an infiltration experiment[END_REF][START_REF] Rizzo | Self-potential signals associated with pumping tests experiments[END_REF][START_REF] Suski | Monitoring of an infiltration experiment using the self-potential method[END_REF][START_REF] Maineult | Self-potential signals induced by periodic pumping tests[END_REF][START_REF] Linde | Self-potential investigations of a gravel bar in a restored river corridor[END_REF] and none of these were instrumented to measure the vertical distribution of the SP signal.

The electro-chemical contribution to the self-potential signal can have two different origins: oxido-reductive (redox) and electro-diffusion (diff) processes. Even if strong redox potential contrasts exist in the near surface, redox phenomena only contribute to the SP signals when an electronic conductor exists that connects two regions with different oxidation potential (e.g., metallic bodies, certain kinds of bacteria). If this is the case, the redox contribution is typically the larger contribution to the SP signal (e.g., [START_REF] Naudet | Groundwater redox conditions and conductivity in a contaminant plume from geoelectrical investigations[END_REF]Linde and Revil, 2007). If not, the corresponding contribution is null (e.g., Hubbard et al., 2011). The electrodiffusive contribution occurs in the presence of concentration gradients and is linked to the differential diffusion of ions with different mobilities. Electro-diffusive phenomena have been extensively studied in saturated porous media (e.g., [START_REF] Maineult | Electrical response of flow, diffusion, and advection in a laboratory sand box[END_REF][START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF][START_REF] Maineult | Influence of the mineralogical composition on the self-potential response to advection of KCl concentration fronts through sand[END_REF][START_REF] Revil | Characterization of transport properties of argillaceous sediments: Application to the Callovo-Oxfordian argillite[END_REF][START_REF] Straface | Estimation of longitudinal dispersivity in a porous medium using self-potential signals[END_REF], but only few works concern partially saturated conditions [START_REF] Revil | Diffusion of ions in unsaturated porous materials[END_REF][START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF]. In the past, the electrodiffusive contribution has often been ignored, for example, during SP monitoring of saline tracer tests (e.g., [START_REF] Bolève | Localization and quantification of leakages in dams using time-lapse self-potential measurements associated with salt tracer injection[END_REF].

We present the first results of a long-term monitoring program designed to investigate the role of SP data for predictive in situ estimation of vertical water flux. The HOBE agricultural test site in Voulund (Denmark) was chosen as the vadose zone is extensive, flow and transport processes can be assumed to be mainly vertical, and it is well instrumented with meteorological, hydrological and geophysical tools and sensors [START_REF] Jensen | HOBE: A hydrological observatory[END_REF]. Vertically distributed non-polarizable electrodes that were installed at the site were monitored for more than 2 years. We first obtain a numerical model of the test site that is used to simulate water fluxes, ionic transport and SP signals. We then compare the predictions for different competing models of SP signal generation and place particular focus on the signal contributions (i.e., electrokinetic and electro-diffusive) during a saline tracer test.

Theoretical framework

Below, we present the theory used to describe water flow, transport and SP signal generation under the assumption of vertical flow and transport only.

Flow and transport

The vadose zone is the region comprised between the land surface and the water table. The water saturation, S w , is defined as the ratio between the water and pore volumes: S w = θ w φ , where θ w is the volumetric water content (m 3 m -3 ) and φ the medium porosity (m 3 m -3 ). The effective water saturation is defined as:

S e = θ w -θ w r φ -θ w r , (1)
where θ w r is the residual water content (m 3 m -3 ). The water retention function relates the effective saturation of the medium to its matric potential, h (m). In this work, we use the van Genuchten (1980) model:

S e = 1+ α VG h ( ) n VG ⎡ ⎣ ⎤ ⎦ -m VG , (2)
where α VG (m -1 ) is proportional to the inverse of the air-entry pressure, while n VG and

m VG = 1-1 n VG ( ) are curve shape parameters.
Water fluxes are described by Richards' equation and the van Genuchten-Mualem model [START_REF] Van Genuchten | A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[END_REF] is used for the relative permeability function, k w rel ,

k w rel (S e ) = S e 1-1-S e 1/m VG ( ) m VG ⎡ ⎣ ⎤ ⎦ 2 . ( 3
)
The hydraulic conductivity as a function of saturation is then given as,

K w (S e ) = k w rel (S e )K w sat , (4)
with (m s -1 ) the saturated hydraulic conductivity. The vertical water flux u (m s -1 ) is described by,

u = -K w (S e ) ∂h ∂z + 1 ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ .
(5)

The pore water is an electrolyte containing N ionic species j with a concentration C j (mol L - 1 ). Transport under partial saturation for each species is driven by the water flux through the 1D conservation equation:

∂ θ w C j ( ) ∂t + ∂ ∂z -α z u θ + D j eff ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ ∂C j ∂z + uC j ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ = 0 , (6)
where α z (m) is the dispersivity of the medium along the z-axis, and D j eff (m 2 s -1 ) is the ionic diffusion coefficient of the j th ionic species in the porous medium.

Self-potential generation

The SP response of a given source current density J S (A m -2 ) can be described by two equations [START_REF] Sill | Self-Potential Modeling from Primary Flows[END_REF],

J = σ E + J S , (7)
∇ ⋅ J = 0 , (8)
where J (A m -2 ) is the total current density, σ (S m -1 ) the bulk electrical conductivity of the medium, E = -∇ϕ (V m -1 ) the electrical field, and ϕ (V) the electrical potential. In absence of external source currents (i.e. no current injection in the medium), Eqs. ( 7) and ( 8) can be combined to obtain:

K w sat ∇ ⋅ σ∇ϕ ( ) = ∇ ⋅ J S .
(9)

The measured SP response at the i th electrode is the potential difference with respect to the reference electrode:

SP = ϕ i -ϕ ref .
In this work, we only consider electrokinetic (superscript EK) and electro-diffusive (superscript diff) contributions. Considering only 1D vertical variations, these sources can be summed to obtain the total source current density:

J S = J S EK + J S diff .
The electrokinetic source ( J S EK ) is directly related to the water flux. This source current density is traditionally defined with respect to the hydraulic gradient through the streaming potential coupling coefficient [START_REF] Helmholtz | Studien über electrische Grenzschichten[END_REF], C EK (V m -1 ):

J S EK = σ C EK ∂h ∂z . ( 10
)
No model has so far been able to predict its dependence on water saturation, C EK (S w ) , for all published data (e.g., [START_REF] Guichet | Streaming potential of a sand column in partial saturation conditions[END_REF][START_REF] Revil | Streaming potentials in two-phase flow conditions[END_REF][START_REF] Allègre | Streaming potential dependence on water-content in Fontainebleau sand[END_REF][START_REF] Vinogradov | Multiphase streaming potential in sandstones saturated with gas/brine and oil/brine during drainage and imbibition[END_REF]).

An alternative approach to C EK proposed by Revil and Leroy (2004) for water saturated conditions was later extended to partial saturation by Linde et al. (2007a). It is based on the excess charge density, Q v (C m -3 ), in the pore water that helps to counterbalance the electrical charges at mineral surfaces (i.e. the charge fraction located in the so-called Gouy-Chapman diffuse layer, Fig. 1a). This excess charge distribution follows a Boltzmann distribution with properties that depend mainly on the pore water salinity and the electrical potential at the mineral surface [START_REF] Hunter | Zeta Potential in Colloid Science: Principles and Applications[END_REF]. When the salinity increases, the thickness of the diffuse layer shrinks (Fig. 1b). When the water flows in the pores (Fig. 1c), the excess charge is dragged in the medium and generates a corresponding current source density. From a variable change in Eq. ( 10), the electrokinetic source current density can be defined by:

J S EK = Q v eff u , (11)
where Q v eff is the excess charge (in C m -3 ) that is effectively dragged in the medium by the water flux u. This concept offers also an alternative definition of the streaming potential coupling coefficient [START_REF] Revil | Constitutive equations for ionic transport in porous shales[END_REF]:

C EK = - Q v eff σ K w ρ w g , ( 12
)
where ρ w (kg m 3 ) is the water density and g the gravity acceleration (9.81 m s -2 ). Titov et al. (2002) presented laboratory data that strongly suggest an inverse relationship between the effective excess charge, , and permeability, k (m 2 ). [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF] proposed the following empirical expression for this type of relationship

log 10 Q v eff,sat ( ) = -0.82 log 10 k ( ) -9.23 . (13)
This relationship has been successfully applied in several studies (e.g., [START_REF] Bolève | Preferential fluid flow pathways in embankment dams imaged by self-potential tomography[END_REF][START_REF] Revil | Coupled hydromechanical and electromagnetic disturbances in unsaturated porous materials[END_REF], even if the influence of the pore water chemistry on Q v eff,sat is neglected. The effective excess charge depends on the salt concentration in the water phase C w (mol m -3 ) as indicated in Fig. 1b (e.g., Hunter, 1981; Revil et al., 1999). This relationship has been observed experimentally for measurements of streaming potential coupling coefficients at different pore water conductivity, σ w (S m -1 ) (i.e., salinity), by Pengra et al. (1999, their Fig. 9). Linde et al. (2007b) propose an empirical relationship to correct for the salinity effect at saturation that is based on literature data for various media:

log 10 C EK sat ( ) = a + b log 10 σ w ( )+ c log 10 σ w ( ) ⎡ ⎣ ⎤ ⎦ 2 ( 14
)
where a = -0.895, b = -1.319, and c = -0.1227 are best fit parameters. Based on a volume averaging approach, Linde et al. (2007a) consider that as the volume of water diminishes in the medium (i.e., the saturation decreases), and given that the amount of surface charge stays constant, the excess charge density Q v eff (S w ) should increase in the pore water. Linde et al. (2007a) propose the following scaling relationship:

Q v eff,sat
Q v eff (S w ) = Q v eff,sat S w . ( 15
)
This relationship provides satisfactory predictions when used to interpret laboratory SP data from drainage and imbibition experiments in rather homogeneous media (e.g., Linde et al., 2007a;[START_REF] Mboh | Coupled hydrogeophysical inversion of streaming potential signals for unsaturated soil hydraulic properties[END_REF][START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF]. Nevertheless, this model is limited when applied to heterogeneous media, especially at low saturations. Indeed, the volume-averaging is only strictly valid when all pores have the same size (see [START_REF] Linde | Self-potential investigations of a gravel bar in a restored river corridor[END_REF][START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF]. [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF] developed two flux averaging approaches that conceptualize the porous medium as a bundle of capillaries (see also, [START_REF] Linde | Comment on "Characterization of multiphase electrokinetic coupling using a bundle of capillary tubes model[END_REF][START_REF] Jackson | Multiphase electrokinetic coupling: Insights into the impact of fluid and charge distribution at the pore scale from a bundle of capillary tubes model[END_REF][START_REF] Jackson | On the validity of the "Thin" and "Thick" double-layer assumptions when calculating streaming currents in porous media[END_REF]. [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF] upscale the electrokinetic properties of a given capillary with radius R that contains a pore water with a given solute concentration, C w , (i.e., the effective excess charge in the capillary Q v eff,R (R,C w ) ) to a representative elementary volume of a porous medium characterized by a saturated capillary size distribution that varies with saturation. This distribution f D can either be derived from the water retention function (referred to as the WR approach in the following) or from the relative permeability function (referred to as the RP approach in the following). The effective excess charge is obtained by integrating the distribution of the pore water flux v R (R) (m s -1 ) and the distribution of excess charge density Q v eff,R (S w ,C w ) (C m -3 ) within the capillaries:

Q v eff (S w ,C w ) = Q v eff,R (R,C w )v R (R) f D (R)dR R min R Sw ∫ v R (R) f D (R)dR R min R Sw ∫ . ( 16
)
Flux averaging based on the WR or the RP approach (

f D WR or f D RP ) yield two different effective excess charge functions (Q v eff,WR or Q v eff,RP
). Previous studies [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF][START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF] indicate that the RP approach reproduce the observed amplitudes well, but is less accurate in describing the relative changes with saturation. The WR approach tends to underestimate the observed amplitudes, but reproduces relative changes with saturation rather well. In analogy with the relative permeability function (i.e., Eq. ( 4)), we decompose the effective charge function into a value at saturation, Q v eff,sat , and a relative saturation-dependent function, Q v eff,rel (S w ) :

Q v eff (S w ) = Q v eff,rel (S w )Q v eff,sat . ( 17
)
The value of Q v eff,sat can be obtained by measurements, by using an empirical law (e.g. Eqs.

(13) or ( 14)), or by inversion. The Q v eff,rel (S w ) is here derived using either the RP or the WR approach (see Eq. ( 16)).

The electro-diffusive source current density, J S diff , is generated by ionic charge separation between anions and cations with different mobilities β j . It can be described by the microscopic Hittorf number which represents the fraction of the total current transported by a given ionic species, j:

t j H = β j β ι ι=1 Q ∑ . ( 18
)
Charged porous media can act as semi-permeable membranes that enhance or decrease this effect and this necessitates a macroscopic Hittorf number T j H . [START_REF] Revil | Diffusion of ions in unsaturated porous materials[END_REF] describe how T j H varies with water saturation and proposed an application to clay-rock. The electro-diffusive source current density is given by,

J S diff = -k B T T j H (S w ) q j σ (S w ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 1 C j ∂C j ∂z j=1 Q ∑ , ( 19
)
where k B = 1.38065 × 10 -23 J K -1 is the Boltzman constant, T (K) is the temperature, q j = ±Z j e 0 (C) is the electrical charge of the considered ions, with Z j its valence and e 0 = 1.3806 × 10 -19 C the elementary charge. The resulting electrical potential is referred to as electro-diffusive or junction potential (e.g., [START_REF] Maineult | Detection of advected concentration and pH fronts from self-potential measurements[END_REF][START_REF] Maineult | Influence of the mineralogical composition on the self-potential response to advection of KCl concentration fronts through sand[END_REF][START_REF] Jouniaux | Review of self-potential methods in hydrogeophysics[END_REF], or membrane potential if the effect of the electrical double layer cannot be neglected (e.g., [START_REF] Revil | Characterization of transport properties of argillaceous sediments: Application to the Callovo-Oxfordian argillite[END_REF]. [START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF] showed that a saturation-independent microscopic Hittorf number t j H (i.e., T j H (S w ) = t j H in Eq. ( 19)) could explain experimental laboratory data.

From Eqs. ( 9), ( 11), and ( 19), the electrical problem can be simplified and re-written as:

∂ϕ ∂z = 1 σ (S w ) Q v eff (S w )u -k B T T j H (S w ) q j σ (S w ) ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 1 C j ∂C j ∂z j=1 Q ∑ ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ . ( 20
)
The superposition principle allows decomposing the resulting SP signal into electrokinetic and electro-diffusive contributions SP = SP EK + SP diff (e.g., [START_REF] Jougnot | Self-Potentials in partially saturated media: the importance of explicit modeling of electrode Effects[END_REF].

Electrical and dielectric petrophysical relationships

In hydrology, dielectric permittivity, ε (F m -1 ), and electrical conductivity, σ (S m -1 ), are often measured to infer soil water saturation and salinity (e.g., [START_REF] Huisman | Measuring soil water content with ground penetrating radar[END_REF][START_REF] Friedman | Soil properties influencing apparent electrical conductivity: A review[END_REF][START_REF] Vereecken | On the value of soil moisture measurements in vadose zone hydrology: A review[END_REF][START_REF] Laloy | Electrical resistivity in a loamy soil: Identification of the appropriate pedo-electrical model[END_REF]. Many petrophysical relationships exist to relate dielectric permittivity and water content (see [START_REF] Huisman | Measuring soil water content with ground penetrating radar[END_REF] for a review).

These relationships are often expressed in term of relative permittivity: ε r = ε ε 0 with ε 0 = 8.854 × 10 -12 F m -1 being the dielectric permittivity of vacuum. Among the existing relationships, one can distinguish between empirical relationships (e.g., [START_REF] Topp | Electromagnetic determination of soil water content: Measurements in coaxial transmission lines[END_REF] and those based on mixing laws or up-scaling procedures (e.g., the Complex Refractive Index Model: [START_REF] Dobson | Microwave dielectric behavior of wet soil-Part II: Dielectric mixing models[END_REF][START_REF] Roth | Calibration of time domain reflectometry for water content measurement using a composite dielectric approach[END_REF]. These more theoretical relationships are of interest as they can be used to explicitly account for the dielectric permittivities of the different components, the medium porosity, and its pore space geometry. In addition, temperature effects on the relative permittivity of pore water can be taken into account by the following empirical relationship [START_REF] Weast | CRC handbook of chemistry and physics[END_REF]: 21) [START_REF] Linde | Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data[END_REF] extended the volume averaging approach of [START_REF] Pride | Governing equations for the coupled electromagnetics and accoustics of porous media[END_REF] to describe the relative permittivity in partially saturated conditions:

ε w (T [°C]) = 78.54 1-4.579 × 10 -3 T -25 ( ) + 1.19 × 10 -5 T -25 ( ) 2 -2.8 × 10 -8 T -25 ( ) 3 ⎡ ⎣ ⎤ ⎦ .(
ε r = φ m S w n ε w + φ -m -1 ( ) ε s + 1-S w n ( ) ε a ⎡ ⎣ ⎤ ⎦ , (22)
where ε w , ε s , and ε a are the relative permittivities of water, grain and air, respectively. The petrophysical parameters m and n are the cementation and the saturation index defined by [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]. They describe the pore space and the fluid phase geometry (i.e., tortuosity and constrictivity), respectively [START_REF] Revil | Electrokinetic coupling in unsaturated porous media[END_REF]. [START_REF] Linde | Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data[END_REF] also derived a petrophysical relationship to predict electrical conductivity:

σ = φ m S w n σ w + φ -m -1 ( ) σ s ⎡ ⎣ ⎤ ⎦ , (23)
where σ w and σ s are the pore-water and the grain surface conductivities (S m -1 ). This petrophysical relationship has been shown to reproduce variations of electrical conductivity with saturation (e.g., [START_REF] Breede | Joint measurement setup for determining spectral induced polarization and soil hydraulic properties[END_REF][START_REF] Laloy | Electrical resistivity in a loamy soil: Identification of the appropriate pedo-electrical model[END_REF]. When the grain surface conductivity can be neglected (e.g., typically in absence of clay), Eq. ( 23) can be simplified [START_REF] Archie | The electrical resistivity log as an aid in determining some reservoir characteristics[END_REF]: σ = φ m S w n σ w .

Material and methods

In this section, we describe the experimental set-up at the Voulund agricultural field site (Denmark), the strategy used for model calibration, and the numerical implementation of our modeling framework described in section 2.

Voulund test site at HOBE and set-up

The Danish hydrological observatory (HOBE) is focused on the Skjern River catchment (Jensen and Illangasekare, 2011) (Fig. 2a). In this catchment, the Voulund field site was chosen for detailed investigations of surface-aquifer processes in an agricultural environment.

This area that is the focus of the present study is characterized by a relatively flat surface topography. The field site is equipped for meteorological, hydrogeological and geophysical monitoring (Fig. 2b andc). A full description of the observatory and the Voulund field site are available online (www.hobe.dk).

At Voulund, the geology is characterized as a fairly homogeneous sandy soil. The water table is monitored on site and it is located between 5.5 and 6.5 m depth. The soils comprising the vadose zone were characterized by five drilling campaigns between September 2011 and April 2012. 6 m long drill cores were extracted from the site and cut into 7.7 cm long samples. The sediment samples were weighed both after collection, and after 48 hours in an oven to determine the volumetric water content, total porosity and saturation. Additionally, 71 samples from one of the five cores were analyzed for grain size distribution. Over the 6 first meters, the soil is mainly composed by more than 90 % of very fine to very coarse sand.

Although, between 5 and 10 % of silt and clay can be found in layer 1 and 2, and less than 5 % in layer 4 and 5. Details about the drilling campaigns and core sample analysis can be found in [START_REF] Uglebjerg | Assessment of groundwater recharge using stable isotopes and chloride as tracers[END_REF]. Furthermore, laboratory measurements were performed on 100 cm 3 soil cores extracted in the near vicinity of the field site. Retention characteristics at 6 depths down to 2.05 m and the soil hydraulic conductivities at 3 depths down to 0.8 m were determined [START_REF] Vasquez | Profile soil water content measurements for estimation of groundwater recharge in different land uses[END_REF]. In general, the sediments constituting the soil can be divided into 7 different layers (Fig. 2c and Table 1).

Table 1: Gaussian prior for the soil properties of each layer in the MCMC inversion: the mean value is expressed on the first line, while the standard deviation is given in parenthesis.

Layer Depth In the present study, meteorological and hydrological data were used as input to characterize the soil and its state evolution through a parameter inversion procedure. Due to the lateral homogeneity of the soil and the absence of significant topography, hydrological processes are expected to be primarily vertical. We used precipitations data from a ground level rain gauge and potential evapotranspiration calculated from solar radiation, wind speed, air temperature, and relative humidity using the Penman-Monteith equation assuming a 0.15 m vegetation cover [START_REF] Monteith | Evaporation and environment[END_REF]. These data were available as hourly values since 2009.

(-) (-) (m -1 ) a (-) (m s -1 ) a (-) b (-) b 1 0 -0.
Monitoring was performed using 17 5TE and 6 MPS-1 Decagon sensors (www.decagon.com) down to 3 m depth. The 5TE sensors measure the relative dielectric permittivity ( ε r ), the temperature (T in °C) and the electrical conductivity (σ in µS cm -1 ) every 0.25 m from 0 to 3 m, while the MPS-1 sensors measure the matric potential (h in m) at 0.25, 0.50, 0.75, 1.00, 2.00, and 3.00 m depth. These measurements were stored in a Campbell datalogger every minutes since July 2011.

The field site was also instrumented to conduct cross-borehole geophysical monitoring with both electrical resistivity tomography (ERT) and ground penetrating radar (GPR) using a total of 9 boreholes of 6 m depth with PVC tubes installed in a square (Fig. 2b). In the current study, we consider the ERT measurements conducted in a five-borehole configuration. Four boreholes form a 5 m large square with an additional borehole at its center. 24 electrodes were installed in each borehole with a spacing of 0.25 m from 0.25 to 6.00 m depth. In addition, electrodes were placed close to the surface (0.25 m depth) along the diagonals of the square.

Each electrode consisted of a 5 cm wide stainless steel mesh, which was wrapped around a PVC tube. After lowering the PVC tubes into their correct location, the drill holes were backfilled with oven-dry sand from the bottom up while gently vibrating the tubes. This was done to avoid cavities along the PVC tubes and thereby ensuring a good electrical contact between the surrounding formation and the electrodes. ERT measurements were acquired on a daily to weekly basis for approximately one year. The measurement configuration details for injection and potential measurement electrodes are described by Haarder et al. (Accepted).

The 3D ERT data were inverted to obtain a 3D distribution of resistivity using the R3t inversion code [START_REF] Binley | R3t version 1.8. Software and userguide[END_REF]. The area within the five boreholes was discretized using an unstructured tetrahedral mesh with a typical length scale of 0.125 m resulting in 614.698 elements. To obtain a 1D resistivity distribution as a function of depth the median and the first and third quartiles of the resistivity at a given depth was used in order to minimize the influence of inversion artefacts.

The SP monitoring was conducted using 15 non-polarizable Pb-PbCl 2 electrodes [START_REF] Petiau | Second generation of lead-lead chloride electrodes for geophysical applications[END_REF] located at 0.25 to 3.20 m depth with the reference electrode at 7.30 m depth (Fig. 2c).

Two electrodes were installed at every depth level in the shallow part (0.25 m, 0.50 m, 0.75 m, 1.00 m, and 1.45 m), while only one electrode was installed at every depth level in the deeper part (1.90 m, 2.50 m, 3.10 m, and 7.30 m). The deepest electrode was chosen as the reference electrode due to its position below the water table. The SP electrodes were installed on July 20 th 2011 and the monitoring is running continuously since then with measurements stored every 5 minutes.

A saline tracer test was conducted in September 2011 to characterize water infiltration processes in situ and to test the ability of geophysical monitoring to quantify groundwater recharge (Haarder et al., Accepted). On the 14 th , the equivalent of 3 mm of brine ( σ w Tr = 230 mS cm -1 ) was infiltrated within 2 hrs across a 142 m 2 area around the borehole square (Fig. 2b). This large tracer injection zone provides an injection that can be approximated as 1D in the vicinity of the geophysical sensors.

Hydrological model calibration

In a first step, we calibrate a hydrological model using the dielectric permittivity and matric potential measurements. Figure 3a describes how this was done using Markov chain Monte Carlo (MCMC) simulation. The prior distributions of the different parameters were chosen based on core sample analysis (Table 1). The hydrological parameters consist of porosities ( φ), saturated hydraulic conductivities ( K w sat ), residual water content (θ w r ), and van Genuchten parameters (α VG and n VG ) that were assigned individually to each of the 7 layers. The petrophysical parameters depend on the relationship used to link dielectric permittivity ( ε sim )

and water content (θ w sim ). Herein, we used the petrophysical relationship by [START_REF] Linde | Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data[END_REF] see Eq. ( 22)). We included additional petrophysical parameters (m and n) as priors in the MCMC procedure. It was assumed that these petrophysical parameters were the same for all soil layers and that surface conductivity could be neglected. The hydrological forward problem (flow and transport) is solved using the Hydrus1D software [START_REF] Simunek | The HYDRUS code for simulating the one-dimensional movement of water, heat, and multiple solutes in variably-saturated media[END_REF] with certain transport properties (α z = 0.25 m, and D NaCl eff = 10 -9 m 2 s -1 ) assumed known (Haarder et al., Accepted). The 7 layers (see Fig. 2c)

were discretized in 300 cells with appropriate boundary conditions assigned based on meteorological data (precipitation and potential evapotranspiration) and water table depth measured on site. In order to limit the influence of the initial conditions, we conducted simulations over 500 days starting on August 1 st 2010; in the following, we consider the day that the tracer test was conducted to be Day 0. The simulation period is discretized in daily time steps (one meteorological input per day) with a refinement to hourly periods after Day -4 to better constrain the time evolution during the tracer test.

The MCMC inversion procedure uses the measured relative permittivity ( ε meas ) and matric potential ( h meas ) as observables. The DREAM (ZS) algorithm is used with three chains and standard values of algorithmic variables; see [START_REF] Laloy | High-dimensional posterior exploration of hydrologic models using multiple-try DREAM (ZS) and high-performance computing[END_REF] for details. In the inversion, we use a standard log-likelihood function that is proportional to:

SSR = 1 2 ε meas -ε sim σ ε std ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 2 ∑ + 1 2 h meas -h sim σ h std ⎛ ⎝ ⎜ ⎞ ⎠ ⎟ 2 ∑ , ( 24
)
where the subscripts meas and sim indicate measured and simulated responses. Standard deviations of σ ε std = 0.5 and σ h std = 0.25 m were assigned for the relative permittivity and the matric potential. The rather large value for σ h std is due to the rather poor quality of the MPS-1 data. This MCMC results represent a full posterior probability density function of the hydrodynamic parameters that explain the hydrologic measurements ( ε meas and h meas ), but we only consider the best fit model in the following.

Hydrogeophysical simulation

The theoretical framework presented in section 2 was implemented to solve the coupled hydro-electrical problem numerically (Fig. 3b). The hydrologic forward model is solved using the best-fit hydrodynamic parameters (see section 3.2) to obtain the vertical distribution of the water flux (u in m s -1 ), the solute concentration ( C w in mol L -1 ), and the water saturation ( S w )

as a function of time. The pore-water conductivity was computed from the solute concentration σ w (C w ) using the Sen and Goode (1992) relationship. It yields a soil electrical conductivity that depends on the water saturation and the solute concentration: σ (S w ,C w ) based on Eq. ( 23).

The effective excess charge functions Q v eff (S w ,C w ) for each geological layer were calculated using the volume averaging model by Linde et al. (2007a), as well as the RP and WR approaches by [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF] using the hydrodynamic parameters obtained from the inversion procedure (see section 3.2). Q v eff is then calculated for each cell and time step based on the S w and C w distribution at the considered time. The electrokinetic contribution to the source current density J S EK is given by Eq. ( 11). The electro-diffusive contribution, J S diff , is calculated from the solute concentration gradients and the medium's electrical conductivity (Eq. ( 19)).

Based on the distribution of electrical conductivity and the total source current density:

, the electrical problem is solved using a modified version of MaFloT (see maflot.com, [START_REF] Künze | An adaptive multiscale method for density-driven instabilities[END_REF]. This yields the electrical potential vertical distribution, thus the SP signal from Eq. ( 20), at the different times (i.e., SP(z,t) ).

J S = J S EK + J S diff

Results and interpretations

Hydrological calibration results

The MCMC inversion procedure described in section 3.2 was applied for data covering a period of 70 days that included the tracer test injection and the monitoring (from Day -2 to Day +68). We initially considered three different petrophysical relationships in the inversion:

the Topp equation, the CRIM, and the one proposed by Linde et al. (2006, Eq. (22)). The latter was retained for further analysis as it can describe both the dielectric permittivity (Eq. ( 22)) and the electrical conductivity (Eq. ( 23)) with the same medium parameters: m and n.

Figure 4 shows the best-fitting predictions of relative dielectric permittivity at different depths (the average RMSE is 0.55). Large ε meas peaks could not be reproduced as they would correspond to unphysical water contents (i.e., largely superior to porosity). This behavior is attributed to unaccounted salinity effects, similar to those studied by [START_REF] Rosenbaum | Correction of Temperature and Electrical Conductivity Effects on Dielectric Permittivity Measurements with ECH O Sensors[END_REF].

Unfortunately, we could not use their empirical correction models as our permittivities were much lower and the salinities were sometimes higher than the ranges considered by [START_REF] Rosenbaum | Correction of Temperature and Electrical Conductivity Effects on Dielectric Permittivity Measurements with ECH O Sensors[END_REF]. The average RMSE for matric potential is 0.37 m. This poor RMSE is likely due to a poor calibration of the sensors to the studied soil. The model parameters with the best fitting predictions are presented in Table 2. The hydrodynamic parameters are consistent with the laboratory characterization of core samples [START_REF] Vasquez | Profile soil water content measurements for estimation of groundwater recharge in different land uses[END_REF]. The petrophysical parameters, m = 1.38 and n = 1.57, are also in good agreement with literature data for sand and sandy soils (e.g., [START_REF] Friedman | Soil properties influencing apparent electrical conductivity: A review[END_REF]Linde et al., 2007a). In the absence of sensors below 3 m, we assume that these petrophysical parameters can also describe the deeper region (between 3 and 7.3 m). ) are used to correct the simulated results using Eq. ( 21). The mean RMSE between the simulated and the measured relative permittivities is RMSE = 0.54.

T meas

Table 2: Best fit estimate of the soil parameters based on MCMC inversion using the Eq. ( 22) proposed by [START_REF] Linde | Improved hydrogeophysical characterization using joint inversion of cross-hole electrical resistance and ground-penetrating radar traveltime data[END_REF].

Layer (-) (-) (cm -1 ) (-) (cm d -1 ) (-) (-) 1

Independent evaluation of the hydrological model

Based on the estimated petrophysical parameters (section 4.1), we calculated the predicted electrical conductivity distribution as a function of depth (σ (z) ) using Eq. ( 23) based on the simulated water saturation ( S w ) and pore water salinity ( C w ). We used the empirical expression by [START_REF] Sen | Influence of temperature on electrical conductivity on shaly sands[END_REF] to transform the water salinity into water conductivity ( σ w ) while accounting for temperature effects. In the following, we will refer to electrical resistivity, the inverse of electrical conductivity: ρ(z) = 1 σ (z) . To evaluate the hydrological model and the petrophysical parameters obtained by the MCMC inversion, we compared our simulation results with selected electrical resistivity profiles obtain by inversion of 3D crossborehole ERT measurements (Haarder et al., Accepted). To do so, we first simulate electrical resistivities based on our Hydrus 1D modeling results ( ρ(z) ). Corresponding electrical tomograms were calculated by forward simulation and subsequent inversion using the R2 code [START_REF] Binley | DC resistivity and induced polarization methods[END_REF]. These results were then compared to the median, the first and the third quartiles of the models obtained from the field inversions. Prior to the tracer injection (Figure 5a), the general trend of the electrical resistivity distribution is well captured by our simulations. This independent evaluation provides some confidence in the hydrological

SP responses in natural conditions

We now turn our attention to the SP data and consider first the raw data (i.e., no filtering or detrending was applied) on Day -2 prior to the tracer injection (Figure 7). The amplitude ranges from SP = -44.6 mV at 0.25 m depth to -11.2 mV at 3.1 m depth. Similar shapes of the SP distribution have been predicted by numerical simulations (e.g., [START_REF] Linde | Self-potential investigations of a gravel bar in a restored river corridor[END_REF][START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF], but never been measured in situ at more than two depths (see [START_REF] Doussan | Variations of self-potential and unsaturated water flow with time in sandy loam and clay loam soils[END_REF]. The observed variations in the SP data at the same depth (maximum difference ~15 mV at 0.25 m depth) are most probably caused by local heterogeneities.

Prior to the tracer injection, the pore water concentration is assumed to be homogeneous ( σ w = 200 µS cm -1 from the first drilling campaign) and the electro-diffusive contribution is null (i.e., J S diff = 0 A m -2 ). In the following, we compare different models that describe

Q v eff = f (S w )
. Unfortunately, we did not have access to samples from the drilling campaign for measuring Q v eff,sat . These samples were taken before the initiation of this work and they were destroyed during the grain size analyses. Instead, we relied on the empirical relationships proposed by [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF] and Linde et al. (2007b) (Fig. 6a) to estimate Q v eff,sat . The model of Linde et al. (2007a) predicts that the effective excess charge increases with the inverse of saturation (Eq. ( 15)). The saturated effective excess charge density was calculated using the empirical relationship by [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF]Eq. (13)) by using the layer permeabilities inferred by the inversion (Table 2). These estimates were subsequently divided by the water saturation at each depth to obtain Q v eff (z) . Then, we solved the electrical problem for SP EK (z) . Note that the empirical relationship proposed by Jardani et al. (2007) (Eq. ( 13))

does not account for any salinity-dependence on Q v eff (see section 2.2 and Fig. 6b). The predicted SP signals (Figure 7) decrease in amplitude with depth and the general shape of the simulated SP signals correspond fairly well to the measurements, but the magnitudes are far too low (1 to 2 orders of magnitude and a corresponding RMSE of 23.7 mV). This result is in agreement with previous field-based studies (e.g., [START_REF] Linde | Self-potential investigations of a gravel bar in a restored river corridor[END_REF]. 13) and ( 15) and the predictions are far too low. The red and blue lines correspond to flux-averaging based on the RP approach ( ) and the WR approach ( ) functions, respectively, scaled with the expected value at saturation. The plain lines and the blue point-dashed line present predictions of the WR approach for a = -0.895 and -0.6 in Eq. ( 14) to obtain the value at saturation, respectively (see Fig. 6). RMSE between the model simulations and the measurements are indicated on the corresponding curves.

Next, we evaluated the RP and WR approaches by [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF]. These formulations allow accounting for the salinity dependence on Q v eff (see section 2.2 and Fig. 1). We compare these results with the empirical relationship by Linde et al. (2007b) (Eq. ( 14) and 6a) that describes the effect of pore water conductivity at saturation. Based on the RP and WR approaches, Figure 6b provides the predicted response of how the effective excess charge at saturation varies with the pore water conductivity; [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF] relationship (Eq. ( 13)) is not salinity-dependent. Figure 6c presents the corresponding variations of the streaming potential coupling coefficients. The streaming potential coupling coefficient predicted by the [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF] relationship overlap well with Eq. ( 14) in a narrow range of pore water conductivities (roughly σ w ∈ 0.01 ; 0.5 [ ] S m -1 ). The salinity dependence of C EK trend is well captured by the RP and WR approaches for σ w > 0.03 S m -1 (Fig. 6c). In the following, we scale the calculated relative effective excess charge functions obtained by the RP and WR models with the prediction at saturation based on Eq. ( 14). obtained by the WR approach using a = -0.6 in Eq. ( 14) as shown in Figure 6a. In the following, we consider only the WR approach scaled with the saturated value from Eq. ( 14)

Q v eff (S w ) f D RP f D WR Fig.
with a = -0.6.

SP responses during the tracer test

The infiltration of the high concentration electrolyte generates Na + and Cl -concentration gradients, which results in differential diffusion (see section 2.2) and a non-negligible electrodiffusive contribution ( J S diff ≠ 0 A m -2 in Eq. ( 20)).

Figure 8a-d shows the raw SP data at the same four times as in Figure 5 together with the predictions based on the WR approach. The influence of the saline tracer injection is evident in the raw data (Fig 8b,c,and d). The SP signal magnitudes diminish, as predicted by theory ( Q v eff diminishes and σ w increases with an increasing salinity); this behavior is clearly seen in the SP model predictions. The raw data also display time-varying gradients between 0 and 1.5 m depth that are likely caused by electro-diffusion. Even if the modeled electro-diffusive contributions have the right shape, we find that the magnitudes are far too low. To better highlight the modeled shape, we have in Figure 8e-h enhanced the simulated electro-diffusive contribution with a factor of 7. This enhancement factor is chosen in the following to obtain magnitudes in accordance with the field data and it has no physical meaning. The simulated water flux and ionic concentrations are displayed in Fig. 9a andb, respectively, and the corresponding simulations of the electrokinetic ( SP EK ) and electro-diffusive ( SP diff ) contributions are displayed in Fig. 9c andd.

The electrokinetic behavior (Fig. 9c) is clearly different before (Day -2) and after (Days +20, +25, and +40) the tracer injection. The saline tracer increases the pore water salinity, and consequently the electrical conductivity. This induces an amplitude diminution of the potential amplitude (Eq. ( 20)), which is enhanced by the resulting diminution of Q v eff (see Figs. 1 and7), and thus also a reduction of J S EK (Eq. ( 11)). This explains why SP EK is almost constant where the tracer concentration is very high (e.g., between 0 and 1.2 m depth on day +20). The electro-diffusive contribution to is calculated from the right term in Eq. ( 20) by using the simulated NaCl concentration ( and from Hydrus 1D transport solution), the water saturation, and the electrical conductivity. In order to determine the macroscopic Hittorf number, we tested two models: from [START_REF] Revil | Diffusion of ions in unsaturated porous materials[END_REF] and the simpler used by Jougnot and Linde (2013) (i.e., the same microscopic Hittorf number as in the pore water). The [START_REF] Revil | Diffusion of ions in unsaturated porous materials[END_REF] model predicts a negative contribution, which is in contradiction with the experimental data (Fig. 8f, g, andh). This might be due to an over-estimation of the influence of the electrical double layer and its effect on the anion mobility that results in anionic exclusion effects. The microscopic Hittorf number predicts a positive contribution, but the predicted signals were (as stated above)

approximately 7 times smaller than the experimental data (Fig. 8b, c, andd). The underlying reasons for this discrepancy is not well understood at the moment. It can partly be attributed

SP diff C Na C Cl T j H (S w ) T j H (S w ) = t j H SP diff SP diff
to differences between our electrical resistivity predictions close to the surface (between 0 and 1.5 m depth) and those obtained by the ERT inversions (see Fig. 5). It is in this region that the concentration gradients are the largest and the electro-diffusive contribution is the most important. We can also not exclude that our calculation of the source term in eq. ( 20) is only partially valid in partially saturated media. This subject has not been addressed in detail in the literature and would be an interesting subject for further studies. We can also not exclude that small-scale heterogeneity, 3-D effects, and preferential flow paths that are not considered in our hydrological model are partly responsible for this discrepancy.

Discussion

In this work, we compared self-potential field data with modeling results in a partially saturated soil both before and during saline tracer infiltration. Prior to the saline tracer injection, we find that the electrokinetic contribution can explain the raw self-potential SP data quite well. This is the first such field demonstration that has been presented in the literature. As the tracer propagates in the medium, a second source type (electro-diffusion) is needed to explain the measurements. For a downward flux, the (Fig. 9c) can only be negative, and a positive signal contribution (such as , Fig. 9d) is needed to explain the change in shape of the vertical self-potential profile when the tracer is present in the soil (between 0 and 2 m depth: Figs. 8f, g, and h).

It is clear that the increase of Q v eff with decreasing saturation predicted by Linde et al. (2007a) is insufficient to explain the observed SP amplitudes (Fig. 6). This finding is consistent with recent studies (e.g., [START_REF] Linde | Self-potential investigations of a gravel bar in a restored river corridor[END_REF][START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF]. The RP and WR approaches proposed by [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF] provide significantly better predictions when scaling the relative excess charge function ( Q v eff,rel (S w ) ) with a value at saturation Q v eff,sat in analogy with the relative permeability function (Eq. ( 4)). This yields simulated results that reproduce the vertical distribution of the electrokinetic contribution. Prior to the tracer injection (Day -2), the SP generation can thus be simulated by only considering the electrokinetic contribution ( from Eq. ( 11)).

To determine a correct Q v eff,sat , we had to account for the influence of salinity upon the excess charge density in the pore water. The widely used approach to determine Q v eff,sat from SP EK SP diff SP EK permeability (Eq. ( 13)) proposed by [START_REF] Jardani | Tomography of the Darcy velocity from self-potential measurements[END_REF] does not account for the salinity effects illustrated in Fig. 1. Figure 6c shows that the saturated streaming potential coupling coefficients obtained from Eq. ( 13) overlap well with the empirical relationship proposed by Linde et al. (2007b; Eq. ( 14)) for typical pore water conductivities (σ w ∈[0.01 ; 0.5] S m -1 ).

This suggests that this empirical relationship is useful for environmental studies under typical salinities, but is of somewhat limited value at low and high salinities. The salinity dependence of C EK trend is well captured by the RP and WR approaches over a wider salinity range when σ w > 0.03 S m -1 .

The SP data are also strongly influenced by the saline tracer, but the present theory is unable to accurately simulate the resulting magnitudes. When the tracer propagates in the medium, the strong concentration gradients make the electro-diffusive contribution quite significant. It

is still an open question how to best model the electro-diffusive contribution under partially saturated conditions. By using a saturation-independent Hittorf number (its value in the pore water), we explain the shape of the SP data after tracer injection (Fig. 8f, g, andh), but the predicted amplitudes appear to be off by a factor 7. This discrepancy is partly attributed to the fact that our predicted electrical conductivities are too high in the tracer-occupied region (Fig. 5b, c, andd). We can also not exclude 3D effects related to the saline tracer plume. Other possibilities relate to the interaction between the salt and the soil matrix or to trapping effects as proposed by [START_REF] Maineult | Influence of the mineralogical composition on the self-potential response to advection of KCl concentration fronts through sand[END_REF] to explain electro-diffusive signal (SP diff ) amplitudes larger than those predicted by theory. Furthermore, only few works have focused on the salinity effect on the Hittorf numbers (e.g., [START_REF] Gulamali | Self-potential anomalies induced by water injection into hydrocarbon reservoirs[END_REF] and the only model to study the saturation effects were proposed by [START_REF] Revil | Diffusion of ions in unsaturated porous materials[END_REF]. This model did not provide results consistent with our data. [START_REF] Jougnot | Derivation of soil-specific streaming potential electrical parameters from hydrodynamic characteristics of partially saturated soils[END_REF] has clearly highlighted the tremendous influence of saturation on streaming potential magnitudes and we suggest that the situation might be similar for the electro-diffusive contribution. This is an important research area to address in future theoretical and laboratory works.

Conclusion

At the Voulund agricultural test site, Denmark, we carried out the first ever field-based monitoring of the vertical distribution of SP signals before and during a saline tracer test. We first derived a hydrological model by MCMC inversion to enable hydrological predictions in response to precipitation and tracer injections. We then propose a hydrogeophysical modeling framework that accounts for both water saturation and salinity variations on the simulated SP signals. This is accomplished by using the concept of an effective excess charge, which varies with water saturation and salinity. The most satisfying modeling results were obtained by the so-called WR approach that conceptualizes the porous media as a bundle of capillaries with a distribution that is inferred from the water retention function. Prior to the tracer injection, we find that the electrokinetic (related to water flow) contribution can explain both the signal amplitude and vertical distribution of the raw self-potential data. This was possible without processing or correcting our field SP data. After the tracer injection, it is clear that an electrodiffusive contribution (related to differential diffusion at concentration gradients) is present.

The predicted shape of the electro-diffusive contribution is in agreement with the field-data, but the magnitudes are too low. We suggest that this is due to inadequacies in the hydrological model and in our petrophysical relationships. This work confirms many theoretical and laboratory findings in showing that the self-potential data is sensitive to water fluxes and concentration gradients. Our theoretical framework is, in principle, able to predict these contributions and show a satisfactory agreement with field data. Nevertheless, accurate predictions of signal magnitudes are complicated by the many subsurface properties (water content, salinity, porosity, pore size distribution, etc.) that affect the data. Future work will focus on strategies to infer long-term infiltration and groundwater recharge at the Voulund agricultural test site using two-years of high-quality monitoring data.

Figure 1 :

 1 Figure 1: Effect of the NaCl electrolyte concentration upon the distribution of excess charge in a capillary with a 40 nm radius: (a) sketch of the electrical double layer, (b) excess charge and (c) pore water velocity distribution. A pore water salinity change from C w = 10 -2 to 10 - 3 mol L -1 result in the Debye lengths increasing from λ D = 3.06 to 9.70 nm. The electrokinetic coupling takes place in the diffuse Gouy-Chapman layer, which has a thickness on the order of two Debye lengths.

Figure 2 :

 2 Figure 2: (a) Location of the Voulund test site in the Skjern river catchment, Denmark (modified from Jensen and Illangasekare, 2011). (b) Overview of sensor positioning and (c) sensor depths. The blue dashed lines in (c) correspond to the maximum and minimum observed water level variation during the study.

Figure 3 :

 3 Figure 3: (a) Hydrological model calibration and (b) hydrogeophysical simulation schemes used in this work.

Figure 4 :

 4 Figure 4: Comparison between the measured (empty circles) and simulated relative dielectric permittivity (filled dots) using the best fitting parameters from the hydrological inversion procedure. The temperatures measured by the 5TE () are used to correct the simulated results using Eq. (21). The mean RMSE between the simulated and the measured relative permittivities is RMSE = 0.54.

  inferred petrophysical parameters. It also indicates that the simulated resistivity magnitudes are correct, which is crucial to accurately simulate SP magnitudes. After the tracer injection (Figure 5b-d), the hydrological simulations indicate a larger decrease of electrical resistivity in the tracer-affected region compared to the ERT results. It is wellknown that ERT provides overly smooth images and it is difficult to assess to which extent the discrepancies in the upper 2 m are related to deficiencies in our hydrological model or to the limited resolution of the ERT results. A certain degree of care is needed when further interpreting results in this depth range.

Figure 5 :

 5 Figure 5: Comparison between inversion results describing the averaged resistivity models with depth. The red dots represent the results obtained from the best fitting hydrological simulation, while the black lines represent the field data inversions from cross-borehole ERT measurements: the solid line is the median of models while the dashed lines are the first and third quartiles of the models. The agreement is good, but the predicted models over-predict the reduction in resistivity due to the tracer injection. The relative RMSE (rRMSE) between the simulated and ERT model (normalized by the ERT model) are indicated for each day.

Figure 6 :

 6 Figure 6: (a) Streaming potential coupling coefficient as a function pore water conductivity in saturated conditions (modified from Linde et al., 2007b). (b) Effective excess charge and (c) streaming potential coupling coefficient as a function of pore water conductivity in saturated conditions relatively to their respective value for = 0.1 S m -1 calculated using the different models in layer 4 and comparison to the models displayed in (a).

Figure 7 :

 7 Figure 7: Comparison between the measured SP signal before the tracer injection and predictions based on different models describing . The black dashed line is based on the volume-averaging approach of Linde et al. (2007a) using Eqs. (13) and (15) and the predictions are far too low. The red and blue lines correspond to flux-averaging based on the RP approach ( ) and the WR approach ( ) functions, respectively, scaled with the

Figure 7

 7 Figure 7 provides a comparison of the SP data with the simulation results based on different

Figure 8 :

 8 Figure 8: Comparison between the measured and simulated SP signal prior (Day -2) and during (Day +20, +25, and +40) the tracer test: the black dots correspond to the measured data and each column to a day. The solid lines in the upper row (a, b, c, and d) correspond to the raw simulated responses (i.e. without enhancement factor), while in the lower row(e, f, g, and h) an enhancement factor of seven for the electro-diffusive contribution was applied. The RMSE between model simulations and measurements are indicated for each subplot.

Figure 9 :

 9 Figure 9: Vertical distribution of the (a) water flux, (b) pore water concentration, and the (c) electrokinetic and (d) electro-diffusive contributions to the total SP signal. Note that the enhancement factor of seven is applied to the electro-diffusive contribution (see Fig. 8f, g, and h).

Acknowledgement

The authors thank the Danish hydrological observatory HOBE for the access to the site, technical help (especially Lars Rasmussen), and full access to the data. The authors used a version of Maflot (maflot.com) that was kindly provided by I. Lunati and R. Künze and modified by M. Rosas-Carbajal. The authors kindly thank the associated editor J.A.

Huismann, S. Ikard, and two anonymous reviewers for their very constructive comments.