Julien Schorsch
email: julien.schorsch@univ-lorraine.fr

Hugues Garnier
email: hugues.garnier@univ-lorraine.fr

Marion Gilson
email: marion.gilson@univ-lorraine.fr

Instrumental variable methods for identifying partial differential equation models of distributed parameter systems

This paper presents instrumental variable methods for identifying partial differential equation models of distributed parameter systems in presence of output measurement noise. Two instrumental variable-based techniques are proposed to handle this continuous-time model identification problem: a basic one using input-only instruments and a more sophisticated refined instrumental variable method. Numerical examples are presented to illustrate and compare the performances of the proposed approaches.

INTRODUCTION

A large variety of natural, industrial, and environmental systems involve phenomena that are continuous functions not only of time, but also of other independent variables, such as space coordinates. Typical examples are transportation phenomena of mass or energy, such as heat transmission and/or exchange, humidity diffusion or concentration distributions. These systems are intrinsically distributed parameter systems (DPS) whose description usually requires the introduction of partial differential equations (PDE). When performing simulations based on a PDE model, it is assumed that all involved parameters are known. However, beforehand these parameters have to be determined. Since they are frequently hardly (or even not at all) a priori known, they have to be fitted from measurements. This amounts to the problem of PDE parameter estimation. There are different ways to identify distributed parameter systems. Three main families of methods can be distinguished [START_REF] Kubrusly | Distributed parameter system identification-A survey[END_REF]):

• methods where optimization techniques are used to directly estimate the PDE model parameters; • methods which aim at approximating the DPS (described by a PDE) to a continuous-time (CT) or discrete-time (DT) lumped parameter system (LPS) (described by ordinary differential or difference equation); • methods which reduce the partial differential equation to an algebraic equation.

Many methods of the three families have been developed.

In [START_REF] Liu | Grey-box Identification of Distributed Parameter Systems[END_REF] an approach is proposed to reduce the PDE models described by two or more independent variables to ordinary differential equations (ODE) models described by only one independent variable. PDE model identification based on continuous space-time separation of the signals is also proposed in [START_REF] Li | Modeling of distributed parameter systems for applications-A synthesized review from time-space separation[END_REF]. The input/output signals are in this case expressed in a basis of polynomial functions such as wavelet (Sadabadi et al. [2008a]), Walsh series [START_REF] Tzafestas | Walsh series approach to lumped and distributed system identification[END_REF]) or Fourier series [START_REF] Sadabadi | Parameter estimation of two-dimensional linear differential systems via Fourier based modulation function[END_REF]). It is also possible to pass through a filtering method involving integrals [START_REF] Sagara | Parameter estimation of distributed-parameter systems under noisy measurements[END_REF]). The problem to estimate the partial derivatives is then got around, the derivatives are no longer approached but integrals are obtained by approximation. In [START_REF] Guo | Identification of partial differential equation models for continuous spatio-temporal dynamical systems[END_REF], an approach which combines implicit Adams integration and an orthogonal least squares is developed.

In [START_REF] Ali | Consistent identification of spatially interconnected systems[END_REF] an identification technique is introduced for discretized PDE models. After discretization of the PDE by finite differences, the parameters are estimated by using techniques developed for discrete-time lumped model identification.

The last decade has witnessed a resurgence of interest devoted to identifying parameters in CT models from DT data. Without relying any longer on analogue computers, the present techniques exploit the power of the digital tools. The developed methods have proven successful in many practical applications and are also available in userfriendly toolboxes for Matlab (see e.g. H. Garnier and L. Wang (Eds.) [2008]). Amongst the available identification approaches for CT input-output models, the interest for instrumental variable (IV) methods has been growing in the last years [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-series Analysis: An Introduction for the Student and Practitioner[END_REF], [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models[END_REF], L. Wang and H. Garnier (Eds.) [2012]). The main reason of this increasing interest is that IV methods offer similar performance as extended least squares (LS) methods or other predictionerror-minimization (PEM) methods and provide consistent results even for an imperfect noise structure which is the case in most practical applications. Amongst the available IV estimators, the Simplified Refined Instrumental Variable for Continuous-time systems, denoted by SRIVC from hereon [START_REF] Young | Refined instrumental variable methods of time-series analysis: Part III, extensions[END_REF], [START_REF] Young | Refined instrumental variable identification of continuous-time hybrid Box-Jenkins models[END_REF]), is, in the authors' experience, the most efficient. It is known to yield asymptotically efficient parameter estimates of CT ordinary differential equation in the presence of white output measurement noise.

The main objective in this paper is to develop a SRIVC version for identifying CT partial differential equation from sampled data.

The paper is organized as follows. Section 2 states the identification problem. In Section 3, two methods based on the simple least squares and a basic instrumental variable estimators are recalled. Both are associated to the statevariable filtering method to handle the partial derivative measurement problem. Section 4 presents the main contribution of the paper, that is the extension of the SRIVC method to identify CT PDE models directly from the sampled data. Section 5 provides numerical simulations in order to illustrate and compare the effectiveness of the proposed methods.

PROBLEM FORMULATION

Consider a two-dimensional linear continuous-time system defined by partial differential equations in the single-input, single-output situation as follows :

nx ix=0 nt it=0 a ix,it p ix x p it t ẙ(x, t) = mx ix=0 mt it=0 b ix,it p ix x p it t u(x, t) (1)
where u(x, t) denotes the input variable observed at space point x and at time t and ẙ(x, t) the noise-free output.

(n x , n t ) are order of the system (assumed to be known), n x ≥ m x , m t and n t ≥ m x , m t and a nx,nt = 1. p x and p t denote the differential operator, respectively

p x = ∂ ∂x , p t = ∂ ∂t (2)
This PDE is subjected to initial and boundary conditions. Equation (1) can also be written in the following compact transfer function (TF) form, ẙ(x, t) = G(p x , p t )u(x, t)

where

G(p x , p t ) = B(p x , p t ) A(p x , p t ) (4) with A(p x , p t ) = nx ix=0 nt it=0 a ix,it p ix x p it t (5) B(p x , p t ) = mx ix=0 mt it=0 b ix,it p ix x p it t (6)
Let h be the spatial sampling distance between two observation points, and T s the time sampling period for the system signals. Denote L as the number of observation points, and N the number of time instants. The sampled signals will be denoted as u(x , t k ) and y(x , t k ) and the output observation equation then takes the form,

S : ẙ(x, t) = G(p x , p t )u(x, t) y(x , t k ) = ẙ(x , t k ) + e(x , t k ) (7)
where ẙ(x , t k ) is the sampled value of the unobserved, noise-free output ẙ(x, t) and e(x , t k ) is a zero-mean, DT

u(x, t) P.D.E. ¨ÿ(x, t) E E m E + e(x , t k ) c y(x , t k )
Fig. 1. Data-generating system white noise source which is assumed to be uncorrelated with the input u(x , t k ) (Fig. 1).

The model is given as follows

G : G(p x , p t , θ) = B(p x , p t , θ) A(p x , p t , θ) = mx ix=0 mt it=0 b ix,it p ix x p it t nx ix=0 nt it=0 a ix,it p ix x p it t (8) where θ = [a 00 • • • a 0nt • • • a nxnt-1 • • • a nx-1nt b 00 • • • b 0mt • • • b mxmt ] T (9)
is the parameter vector that characterises the model (θ ∈ R nx * nt+mx * mt-1 ).

Finally, the identification model can be written in the form

M : ẙ(x, t, θ) = G(p x , p t , θ)u(x, t) y(x , t k , θ) = ẙ(x , t k , θ) + e(x , t k ) ( 10
)
The objective is then to estimate the parameters (9) of the model (10) based on the sampled input and output data

Z L×N = {u(x , t k ), y(x , t k )} L,N =1,k=1 .

BASIC LS AND IV-BASED SVF ESTIMATORS

In comparison with discretised PDE model identification, direct identification of PDE model raises several technical issues. Unlike the discretized PDE model, where only sampled input and output appear, the CT PDE contains partial derivative terms that are required, which are usually not available from measurement. This problem is well-known and various methods have been devised to handle the time-derivatives in the case of ordinary differential equation models [START_REF] Garnier | Continuous-time model identification from sampled data. Implementation issues and performance evaluation[END_REF]). One traditional approach that dates from the days of analogue computer is known as the state variable filtering (SVF) method. The main idea is to pre-filter the PDE model so that the pre-filtered partial derivatives can then be exploited for model parameter estimation. Basic LS and instrumental variable (IV) solutions are described in this section.

The traditional LS-based SVF estimator

Let us first consider the identification problem in the noisefree case. Equation ( 3) can be written in the form

A(p x , p t )ẙ(x, t) = B(p x , p t )u(x, t) (11)
Let F (p x , p t ) be a filter applied to (11),

A(p x , p t )F (p x , p t )ẙ(x, t) = B(p x , p t )F (p x , p t )u(x, t) (12)
The minimum-order SVF filter is typically chosen to have the following form

F (p x , p t ) = 1 (p x + λ x ) nx (p t + λ t ) nt (13)
where λ x and λ t are breakpoint frequencies. This filter is n x -order in space and n t -order in time. Let F jt ix (p x , p t ) for i = 0, 1, • • • , n x and j = 0, 1, • • • , n t be a set of filters defined as

F jt ix (p x , p t ) = p ix x p jt t (p x + λ x ) nx (p t + λ t ) nt (14)
when the initial and boundary conditions are both null, by using the filters defined in ( 14), ( 12) can be rewritten as

nx ix=0 nt it=0 a ix,it F it ix (p x , p t ) ẙ(x, t) = mx ix=0 mt it=0 b ix,it F it ix (p x , p t ) u(x, t) (15) or nx ix=0 nt it=0 a ix,it ẙit ix (x, t) = mx ix=0 mt it=0 b ix,it u it ix (x, t) (16)
where

ẙjt ix (x, t) = F jt ix (p x , p t )ẙ(x, t) u jt ix (x, t) = F jt ix (p x , p t )u(x, t) (17)
At time-instant t = t k and for a space point x = x , considering now the situation where there is an additive noise on the output measurement, equation ( 15) can be rewritten as

y nt nx (x , t k ) = ϕ T (x , t k )θ + (x , t k ) (18) with ϕ T (x , t k ) = -y nt nx-1 (x , t k ) • • • -y nt 0 (x , t k ) -y nt-1 nx (x , t k ) • • • -y 0 0 (x , t k ) u mt mx (x , t k ) • • • u 0 0 (x , t k ) (19)
The LS-based SVF parameter estimates are given by

θLSSV F = L =1 N k=1 ϕ(x , t k )ϕ T (x , t k ) -1 • L =1 N k=1 ϕ(x , t k )y nt nx (x , t k ) ( 20
)
provided that the inverse exists.

Input-only IV-based SVF estimator

Although simple, the LS-based SVF estimates are always biased when the observed data are corrupted by noise. A solution to obtain consistent estimates is to use an instrumental variable method [START_REF] Söderström | Instrumental Variable Methods for System Identification[END_REF], [START_REF] Young | Recursive Estimation and Time-series Analysis: An Introduction for the Student and Practitioner[END_REF]).

The instrumental variable estimator is a classical variant of least squares method. The principle of this estimator is to introduce a new vector such as its components called instruments or instrumental variables denoted here as ζ(x , t k ) are correlated with the regression vector ϕ(x , t k ) but uncorrelated with the noise.

E[ζ(x , t k )ϕ T (x , t k )] be not singular (21a) E[ζ(x , t k ) (x , t k )] = 0 (21b)
where E(.) denotes the mathematical expectation.

The common IV-SVF based parameter estimate is then given by

θIV SV F = L =1 N k=1 ζ(x , t k )ϕ T (x , t k ) -1 • L =1 N k=1 ζ(x , t k )y nt nx (x , t k ) (22)
Different choice for the instruments have been developed. A first and simple solution was presented in [START_REF] Sagara | Parameter identification of distributed-parameter systems in the presence of measurement noise[END_REF]. It involves the use of filtered partial derivatives of the input as instruments. The so-called input-only instrument is defined by

ζ T (x , t k ) = -u nt nx-1 (x , t k ) • • • -u nt 0 (x , t k ) -u nt-1 nx (x , t k ) • • • -u 0 0 (x , t k ) u mt mx (x , t k ) • • • u 0 0 (x , t k ) (23)
Since the instrument is built up from filtered input signals assumed to be uncorrelated with the disturbance, the condition (21b) is satisfied. However the choice of the input type might lead to parameter estimates with quite poor accuracy as illustrated in Section 5.

Comments

(1) The SVF is taken for a large observation time T and a space interval [0, X], the terms related to the initial and boundary conditions decay exponentially provided the system is stable and become negligible quite quick. Thus these terms may be neglected after a time t 0 = k 0 T s and a space point x 0 = 0 h. The estimation algorithm is then applied over [x 0 , X] × [t 0 , T ], where t 0 and x 0 have to be chosen comparable to the settling time of the SVF filters. The number of parameters to be estimated can in this way be reduced substantially and this is surely advantageous with regard to computation efforts and numerical properties.

Fig. 2. Estimation area

(2) An another approach to get rid off the initial and boundary condition problem consists in using a modulating function based method (see e.g. [START_REF] Sadabadi | Parameter estimation of two-dimensional linear differential systems via Fourier based modulation function[END_REF]). A clear advantage of the modulating function methods is that because of their properties, the effects of initial and boundary conditions disappear in the estimation model equation.

(3) The user must provide a priori the breakpoint frequencies λ x and λ t of the SVF (13). Intuitively, this can be chosen in order to emphasize the frequency band of interest and generally, they should be chosen equal to, or larger than, the bandwidth of the system to be identified. The a priori choice can affect the variance of the parameter estimates. (4) A second IV method consists in using an instrumental variable built up from an auxiliary model. The instrument is then given by

ζ T (x , t k ) = - ŷnt nx-1 (x , t k ) • • • - ŷnt 0 (x , t k ) - ŷnt-1 nx (x , t k ) • • • - ŷ0 0 (x , t k ) u mt mx (x , t k ) • • • u 0 0 (x , t k ) ( 24
)
where ŷ(x , t k ) is the output of the auxiliary model obtained from the LS-SVF estimate (20). This method, as the previously presented input-only IV gives unbiased parameter estimates, but not minimum variance. (5) The digital implementation of the various continuous filtering operations can have an important impact on the quality of the parameter estimates. A possibility consists in discretizing the continuous-time filters by the bilinear method as proposed in [START_REF] Sagara | Parameter identification of distributed-parameter systems in the presence of measurement noise[END_REF].

REFINED IV ESTIMATOR

For ordinary differential equations, one particularly successful IV-based stochastic identification method is the iterative Simplified Refined Instrumental Variable (SRIVC) method for continuous-time model identification. This method of adaptive prefiltering is based on a quasi-optimal statistical solution to the problem when the additive noise is white. SRIVC is a logical extension of the more heuristically defined SVF (see e.g. [START_REF] Garnier | Time-domain approaches to continuous-time model identification of dynamical systems from sampled data[END_REF], [START_REF] Young | Recursive Estimation and Time-series Analysis: An Introduction for the Student and Practitioner[END_REF] and the references therein). This technique presents also the advantage of not requiring manual specification of prefilter parameters. In the following, the SRIVC method is extended to handle the case of PDE.

Following the usual prediction error minimization (PEM) approach, a suitable error function (x, t) is given,

(x, t) = y(x, t) - B(p x , p t ) A(p x , p t ) u(x, t). ( 25
)
Minimization of a least squares criterion function in (x, t) measured at the sampling instants provides the basis for the output error estimation methods. However (x, t) can also be rewritten as

(x, t) = 1 A(p x , p t ) (A(p x , p t )y(x, t) -B(p x , p t )u(x, t))(26)
Since the operators commute in this linear case, the filter where the set of filters now takes the form

F (p x , p t ) = 1 A(p x ,
F jt ix (p x , p t ) = p ix x p jt t A(p x , p t ) ( 30
)
The associated estimation model can be written at timeinstant t = t k and at space point x = x in the form

y nt nx (x , t k ) = ϕ T (x , t k )θ + (x , t k ) (31)
where ϕ T (x , t k ) and θ are defined as in ( 19) and ( 9) respectively with F jt ix (p x , p t ) now defined in (30). The proposed SRIVC algorithm for PDE model identification is summarized below. It involves an iterative procedure in which, at each iteration, the auxiliary model used to generate the instruments, as well as the associated pre-filter, are updated, based on the parameter estimates obtained at the previous iteration.

• Step 1 Initialisation

Generate an initial estimate of the TF model parameter vector θ0 using for instance the LSSVF estimator and use this to define the initial prefilter

F (p x , p t , θ0 ) = 1 Â(px,pt, θLSSV F ) , iter = 1. • Step 2
Compute an estimate of the noise-free output ŷ(x, t) by simulating the auxiliary model ŷ(x, t) = B(p x , p t , θiter-1 )

Â(p x , p t , θiter-1 ) u(x, t) (32)
based on the estimated parameters θiter-1

• Step 3 Compute F jt ix (p x , p t , θiter-1 ) = p ix p jt A(p x , p t , θ iter-1 ) (33)
Prefilter the input u(x , t k ), the output y(x , t k ) and auxiliary model output ŷ(x , t k ) by the filter F jt ix

•

Step 4 Form the filtered estimated regressor as

ϕ T (x , t k ) = -y nt nx-1 (x , t k ) • • • -y nt 0 (x , t k ) -y nt-1 nx (x , t k ) • • • -y 0 0 (x , t k ) u mt mx (x , t k ) • • • u 0 0 (x , t k ) (34)
with the instrument

ζ T (x , t k ) = - ŷnt nx-1 (x , t k ) • • • - ŷnt 0 (x , t k ) - ŷnt-1 nx (x , t k ) • • • - ŷ0 0 (x , t k ) u mt mx (x , t k ) • • • u 0 0 (x , t k ) (35) • Step 5 Compute the IV estimate θiter SRIV C = L i=1 N k=1 ζ(x , t k )ϕ T (x , t k ) -1 • L i=1 N k=1 ζ(x , t k )y nt nx (x , t k ) (36)
where θiter SRIV C is the SRIVC estimate at iteration iter based on filtered data. If convergence occurs or the maximum number of iterations is reached then stop, else set iter = iter + 1 and go to Step 2.

Comments

(1) The SRIVC method is known to be asymptotically unbiased and optimal if the noise on the output is white in the case of ODE. These properties are therefore expected in the proposed version PDE and will be illustrated in the next section.

(2) The digital implementation of the analog filtering operations can be achieved, as for the SVF-based estimates, by the bilinear method or by using traditional PDE discretization techniques [START_REF] Claes | Numerical Solution of Partial Differential Equations by the Finite Element Method[END_REF]). Filter stability is then provided while the PDE discretization scheme is stable.

NUMERICAL EXAMPLES

In this section, several simulated examples are considered to illustrate and compare the performances of the proposed approaches. In these examples, the unknown parameters are estimated by three methods, i.e. the least squares based SVF (denoted as LSSVF), the input-only IV-based SVF (denoted as IVSVF) and the Refined IV (denoted as SRIVC). The breakpoint frequencies are taken as λ x = 10 and λ t = 10 (see [START_REF] Sagara | Parameter identification of distributed-parameter systems in the presence of measurement noise[END_REF]).

Example 1

This example comes from [START_REF] Sagara | Parameter identification of distributed-parameter systems in the presence of measurement noise[END_REF]. It is used here to evaluate the effectiveness of the input-only IVSVF method. The same simulation conditions are used.

     1 ∂ 2 ẙ ∂t 2 (x, t) -2 ∂ 2 ẙ ∂x 2 (x, t) = u(x, t) 0 ≤ t ≤ 10 0 ≤ x ≤ 2 y(x , t k ) = ẙ(x , t k ) + e(x , t k ) (37)
The input and the noise-free output are u(x, t) = -4(x 2 + 4)cos(t) ẙ(x, t) = 4x 2 cos(t)

The sampling intervals are taken to be h = 0.25, T s = 0.01 (39)

The input (Fig. 3) and the output are sampled at nine space points and 1000 times instants. The measured output (Fig. 4) is corrupted by a white noise with zero-mean and space-dependent variance σ 2 x such that [START_REF] Sagara | Parameter identification of distributed-parameter systems in the presence of measurement noise[END_REF])

σ x = 0.1x 2 (40)
Monte Carlo (MC) simulation results of 100 runs are presented in Table 1 where the mean and standard deviation of the estimated LSSVF and IVSVF parameters are displayed. It can be seen that the LSSVF estimates are, as expected, biased while the IVSVF method delivers unbiased estimates of the model parameters with reasonable standard deviation in this situation where the input is well-correlated with the output. With the simulation conditions chosen, the input is well correlated with the system output which gives unbiased estimates for the IVSVF. However it must be noted that the choice of the breakpoint frequencies can have a considerable effect on the quality of the parameter estimates. Furthermore as illustrated in the following example, the results deteriorate when the input is less correlated with the output.

Example 2

Let us now consider the following PDE model

     1 ∂ 2 ẙ ∂t 2 (x, t) -0.5 ∂ 2 ẙ ∂x 2 (x, t) = u(x, t) 0 ≤ t ≤ 10 0 ≤ x ≤ 1 y(x , t k ) = ẙ(x , t k ) + e(x , t k ) (41)
The input and the noise-free output are

u(x, t) = 6x(1 -x)t -t 3 ẙ(x, t) = x(1 -x)t 3 (42)
In this example, T s = 0.01s and h = 0.01. The noise is taken as a 2-D zero-mean normally distributed white noise.

The results are shown in Table 2 and were obtained from a MC simulation of 200 runs. In this situation, the input is not well-correlated with the output. The LSSVF estimates are again biased since the output is noisy. Moreover, the IVSVF method leads to poor results in this case. This can be explained from the conditioning of the matrix to be inverted (see eq 21a). Indeed with the chosen input, the condition number of the matrix to be inverted is much poorer (1.7 × 10 17 ) than in example 1 where the condition number was 152. The matrix to be inverted becomes then ill-conditioned which leads to erroneous parameter estimates. Unlike the IVSVF method, the proposed SRIVC method does not suffer from this drawback and gives unbiased parameter estimates.

Example 3

To study further the effects of the input type on the performance of the 3 estimators, another system is also considered.

     1 ∂ 2 ẙ ∂t 2 (x, t) -2 ∂ 2 ẙ ∂x 2 (x, t) = u(x, t) 0 ≤ t ≤ 100 0 ≤ x ≤ 1 y(x , t k ) = ẙ(x , t k ) + e(x , t k ) (43)
The input and the noise are taken as 2-D zero-mean normally distributed white noise and the output is obtained here from the discretization of the partial differential equation by a finite difference method. The data set used for identification involves 100×10000 samples for h = 0.01 and T s = 0.01. Table 3 shows the results of a MC simulation of 200 runs. As in example 2, it can be noticed that both LSSVF and IVSVF estimates are again not acceptable, contrary to the SRIVC estimates which are unbiased with low standard deviation.

CONCLUSION

In this paper, a refined IV-type method has been proposed to directly estimate asymptotically efficient parameters in continuous-time partial derivative equation models from sampled data, when the additive measurement noise is white. One of the main advantages of this method is that it provides partial differential equation models whose parameters can be interpreted immediately in physically meaningful terms.

Fig. 3 .

 3 Fig. 3. Input (in Example 1)

Table 1 .

 1 MC simulation results (example 1)

			a 02	a 20
	Method True value	1	-2
	LSSVF	mean	0.0509 -2.1475
		std	0.0166	0.0085
	IVSVF	mean	1.0535 -1.9964
		std	0.3579	0.0331

Table 2 .

 2 MC simulation results (example 2)

			a 02	a 20
	Method	True value	1	-0.5
	LSSVF	mean	0.5824 -1.1202
		std	0.1648	0.2449
	IVSVF	mean	0.6467 -1.1216
		std	3.6778	6.5423
	SRIVC	mean	1.0418 -0.4983
		std	0.2698	0.2284

Table 3 .

 3 MC simulation results (example 3)

			a 02	a 20
	Method	True value	1	-2
	LSSVF	mean	0.14231 -0.13226
		std	0.0009	0.0007
	IVSVF	mean	1.0422	-0.37014
		std	0.0668	0.0031
	SRIVC	mean	0.99614	-1.9943
		std	0.0152	0.0243


Franck Gilbert
email: gilbert@com.univ-mrs.fr

Stefan Hulth

Niklas Strömberg

Katja Ringdahl

Jean-Christophe Poggiale

Niklas Stro

2-D optical quantification of particle reworking activities in marine surface sediments

Keywords: Amphiura filiformis, Luminophores, Sediment reworking, Two-dimensional optical quantification

Introduction

In bioturbated sediments, particles are subject to various uni-or omnidirectional displacements due to reworking activities of benthic infauna [START_REF] Rhoads | Organism -sediment relations on the muddy sea floor[END_REF][START_REF] Aller | The effects of macrobenthos on chemical properties of marine sediment and overlying water[END_REF][START_REF] Berkenbusch | Factors influencing sediment turnover by the burrowing ghost shrimp Callianassa filholi (Decapoda: Thalassinidea)[END_REF]Franc ß ois et al., 2001). Biologically mediated transport of organic and inorganic particles across redox boundaries may considerably affect rates and pathways of organic matter mineralization during early diagenesis in surface sediments [START_REF] Aller | Bioturbation and remineralization of sedimentary organic matter: effects of redox oscillation[END_REF][START_REF] Gilbert | In situ bioturbation and hydrocarbon fate in an experimental contaminated Mediterranean coastal ecosystem[END_REF][START_REF] Hulthe | Effect of oxygen on degradation rate of refractory and labile organic matter in continental margin sediments[END_REF][START_REF] Sun | Enhanced degradation of algal lipids by benthic macrofaunal activity: effect of Yoldia limatula[END_REF]. Due to the strong influence on sedimentary transformation reactions and transport processes, sediment reworking is a key component in the understanding of overall sediment biogeochemistry, and quantification of biologically induced particle transport is a vital parameter for large-and small-scale diagenetic models [START_REF] Soetaert | A model of early diagenetic processes from the shelf to abyssal depths[END_REF][START_REF] Boudreau | Diagenetic Models and their Implementation[END_REF].

Generally, biological particle reworking is evaluated by diffusion models [START_REF] Berner | Early Diagenesis: A Theoretical Approach[END_REF] and quantified by the associated biodiffusion coefficient (D b ) predicted from the model. This coefficient of sediment reworking is normally verified by experimental manipulations such as additions of fluorescent particles [START_REF] Ge ´rino | Comparison of different tracers and methods used to quantify bioturbation during a spring bloom: 234-thorium, luminophores and chlorophyll-a[END_REF], or direct quantification of the vertical distribution patterns of more or less conservative tracer elements, e.g. 210 Pb, 234 Thorium, Chl. a (e.g., [START_REF] Ge ´rino | Comparison of different tracers and methods used to quantify bioturbation during a spring bloom: 234-thorium, luminophores and chlorophyll-a[END_REF][START_REF] Berg | Effects of bioturbation on solutes and solids in marine sediments[END_REF]. From a natural 3-D displacement of luminescent particles or particle adsorbed tracers in the sediment, the distribution of tracer is projected in horizontal layers with a resulting 1-D repartition. Generally, the sediment has been sectioned in cm-thick layers and the experiment lasts at least a few days.

Recent developments of planar optodes (or plate fluorosensors) have allowed direct quantification of 2-D solute distributions of O 2 and pH [START_REF] Glud | Planar optrodes: a new tool for fine scale measurements of two-dimensional O 2 distribution in benthic communities[END_REF][START_REF] Glud | An in situ instrument for planar O optode measurements at benthic interfaces[END_REF][START_REF] Hulth | A pH plate fluorosensor (optode) for early diagenetic studies of marine sediments[END_REF]. Main principles of these optical plate sensors is to screen an illuminated fluorescent solute specific sensor foil in direct contact with the sediment for high spatial (50 -100 Am) and temporal (less than minute) resolution 2-D imaging of solute distributions patterns. Intensity of fluorescence is directly related to solute concentration.

The overall objective of this study was to develop a fast, nondestructive technique for direct determination of 2-D sediment reworking rates with high temporal and spatial resolution. To fulfill these requirements, we made use of a similar optical setup as for the 2-D pH sensor [START_REF] Hulth | A pH plate fluorosensor (optode) for early diagenetic studies of marine sediments[END_REF]. Quantification of sediment reworking rates with this novel optical approach was applied at a site characterized by intense reworking activity mainly caused by a high density population of Amphiura filiformis (up to more than 3500 ind. m À 2 ; [START_REF] Rosenberg | Benthic marine fauna structured by hydrodynamic processes and food availability[END_REF]. This brittle star switches between suspension and deposit feeding depending on hydrodynamic conditions in the bottom water and organic matter availability [START_REF] Buchanan | A comparative study of some features of the biology of Amphiura filiformis and Amphiura chiajei (Ophiurroidea) considered in relation to their distribution[END_REF], and uses its arms in the sediment to adjust its position and to direct its movement [START_REF] Rosenberg | Density-dependent migration in an Amphiura filiformis (Amphiuridae, Echinodermata) infaunal population[END_REF]. A. filiformis is known to be a highly active sediment reworker [START_REF] Rosenberg | Density-dependent migration in an Amphiura filiformis (Amphiuridae, Echinodermata) infaunal population[END_REF], as recently confirmed also by time-lapse sediment profile imagery [START_REF] Solan | Observation and quantification of in situ animal -sediment relations using timelapse sediment profile imagery (t-SPI)[END_REF].

Materials and methods

Sites and sediment sampling

Stations in the Gullmars Fjord and the Skagerrak area were visited in August 1999 and May 2001 (Fig. 1; Table 1). In 1999, sediment cores (i.d. 10 cm, 30 cm height) were sampled at each of the five stations (S1 -S5) using a multiple corer [START_REF] Barnett | A multiple corer for taking virtually undisturbed samples from shelf, bathyal and abyssal sediments[END_REF]. During the 2001 cruise, only the Deep Trench station in the Skagerrak was studied (S5; Table 1). Sediment subsamples were taken from an Olausson box corer (50 Â 50 cm) using transparent Plexiglas rectangular cores (4 cm side, 30 cm height).

Sediment reworking

In 1999, the sediment cores were incubated at in situ temperature under a water flow system allowing three complete turnovers of water each day in each core. In order to quantify the sediment reworking activity in the sediments, the overlying water was removed and fluorescent luminophores (4 g; 80 -150 Am diameter) were deposited at the sediment surface [START_REF] Ge ´rino | The effects of bioturbation on particle redistribution in Mediterranean coastal sediment. Preliminary results[END_REF]. Then, the water flow system was reactivated. After 9 days of incubation, the sediment cores were sectioned and each layer (1 cm thick) was sieved through a 250-Am mesh to separate macrofauna (preserved in 4% buffered formaldehyde for individual identification) and particles. Sediments were then dried and homogenized. Sediment subsamples were taken for luminophore counting under UV light [START_REF] Mahaut | A luminophore tracer technique for bioturbation studies[END_REF].

The diffusive-like transport of luminophores demonstrated by the tracer repartition profiles in the different sediments allowed us to describe mathematically this bio-mediated transport process as a diffusive one. Then, for each core, a biodiffusion coefficient (D b ) was calculated using the Cochran model (1985).

In 2001, the sediment cores were incubated at in situ temperature under a water flow system. Fluorescent luminophores (2 g; 80 -150 Am diameter) were deposited at the sediment surface. Cores were then precisely positioned at various time intervals in front of the optical system and pictures of core faces were taken. The optical system used in this study was based on the system already used by [START_REF] Hulth | A pH plate fluorosensor (optode) for early diagenetic studies of marine sediments[END_REF] for 2-D quantification of pH. Excitation light (300 W Xe UV/VIS arc lamp; band pass 510 nm excitation filter) was transferred to the sample by a liquid light guide. Light emitted by the luminophores was collected, via an emission Hi Pass 460 nm filter (absorbs every wavelength up to 460 nm), by an 8-bit COHU 4910 integrating monochrome CCD camera (chip resolution 753 Â 582 pixels). The theoretical individual pixel resolution for the COHU camera as set up was 8.5 Â 8.5 Am. NIKON standard single-lens reflex (SLR) lenses were used to optimize spatial resolution and size of the image. Effective image resolution was 75 Â 75 Am at the distance employed between core and camera (20 -40 cm). The data acquisition software InCytIm2y was controlled by a Windows-based PII-400 MHz computer.

Images obtained with InCytIm2y software were thresholded and binarized in order to get ''luminophore pixels'' coded with 1 and ''sediment pixels'' coded with 0. Each image was thus associated to a matrix of 0 and 1. The temporal evolution between two images (e.g. M1 and M2) was then globally quantified by the difference of the matrices corresponding to the different times. This process led to a new matrix (called M) of 0, 1 and À 1 values containing a global information of the movement of luminophores with time (Fig. 2). The value 0 was related to the pixels that did not change. The value 1 was related to the pixels that moved from ''sediment'' to ''luminophore'' state, and inversely for the value À 1.

A reworking coefficient was then constructed to summarize the total reworking information contained in the matrix M (Fig. 2). The sum of all the numbers in the M matrix quantifies the total variation of amount of luminophores. A positive result means that the number of values equal to 1 is larger than the number of values equal to À 1. As the sum was image size dependent, we proceeded then as follows:

-for a fixed number k, all the cells of size k  k pixels contained in the image were considered, -a matrix was associated to each cell in the same way as for the complete image, -the sum of the numbers in each matrix and the mean of these sums were calculated, k was varied from 1 to the total image size.

For each k, a number which represented the mean variation of amount of luminophores in cells of size k was obtained (unit: DQ cm À 1 ). This number increases linearly with respect to k. The slope, given by a linear regression, that quantifies the reworking is sizeindependent. In order to get a rate of reworking quantifying the mean movement of luminophores per surface unit and time (DQ cm À 2 t À 1 ), the slope was divided by the time interval corresponding to the initial and the final images. This rate was named ORC, for optical reworking coefficient (Fig. 2).

Results

Quantification of sediment reworking activities by model predictions on 1-D luminophore distribution patterns revealed comparable biodiffusion coefficients (average D b : 14 cm À 2 year À 1 ) and sediment depth of maximum particle burial (averaging 5.3 cm) for the Fig. 2. Two-dimensional optical luminophore method. Schematic diagram based on a simple theoretical 3 Â 3 matrix illustrating the calculation of the rate of reworking (named ORC, for optical reworking coefficient) corresponding to the mean movement of luminophores per surface unit and time (unit: DQ cm À 2 t À 1 ).

Gullmars Fjord stations (S1 -S3; Table 1). The two Skagerrak stations (S4 and S5) were distinctly different from the Gullmars Fjord stations, and individually separated. Sediment reworking at the Norwegian Trench station (S4) was the lowest measured (D b = 5.4 F 0.6 cm À 2 year À 1 ; n = 3), while calculated coefficient of biodiffusion was the highest measured (D b = 35.5 F 3.7 cm À 2 year À 1 ; n = 3) at the Deep Trench (S5). At the latter station, fluorescent particles were found down to about 7 cm (6.7 F 3.7), while at S4 transport of particles was not significant below a sediment depth of 3 cm (Table 1).

Different sets of images were obtained with the optical system. As presented in Fig. 3, images sequences clearly showed rapid nonlocal transport of fluorescent particles initially deposited at the sediment surface.

Fig. 4 presents six examples, among the total set of 24, of optical reworking coefficients and associated regression coefficient (R 2 ) calculated from Skagerrak Deep Trench station (S5) sediments. Values ranged from 13.95 to 44.53 DQ cm À 2 h À 1 with a mean ( F S.D.; n = 24) value of 27.42 ( F 9.01) DQ cm À 2 h À 1 . Associated regression coefficients were never less than 0.95.

Closer study of the images taken from the four faces of two cores showed that the calculated ORC was variable between the core faces, e.g. being up to 2.5 times higher in one face compared to another (core A; Table 2). However, the mean ORC calculated from results of the four faces of each core were almost similar: 26.85 F 10.45 DQ Fig. 4. A. filiformis (S5) and two-dimensional optical luminophore method. Six examples of calculation of the optical reworking coefficient (ORC; unit: DQ cm À 2 h À 1 ) and associated regression coefficient (R 2 ). cm À 2 h À 1 (F S.D.; n = 4) and 22.63 ( F 6.82) DQ cm À 2 h À 1 for core A and B, respectively.

A long-term image sequence was performed from one face of the third core. Resulting optical reworking coefficients (ORC) were calculated either ''from start'', using the first image of the sequence as reference (Dt between image i and first image), or ''step-by-step'', using two consecutive images (Dt between image i + 1 and image i). The ''from start'' calculation showed a stable ORC (f 27 DQ cm À 2 h À 1 ) with time, excepted a slight decrease down to f 14.87 DQ cm À 2 h À 1 after 38.45 h of incubation (Fig. 5). On the other hand, the ''step-by-step'' calculation demonstrated only a general tendency of the ORC to be relatively more spread with time (from 23.04 to 28.11 DQ cm À 2 h À 1 ). Optical reworking coefficients (ORC) and associated regression coefficient (R 2 ) calculated with two images (Dt: from 8 to 12 h), for each face of two cores. For the method of ORC calculation, see text.

Fig. 5. A. filiformis (S5) and two-dimensional optical luminophore method. Optical reworking coefficient (ORC) vs. time, calculation either ''from start'', using the first image (t = 0) of the sequence as reference (Dt between image i and first image; solid symbols), or ''step-by-step'', using two consecutive images (Dt between image i + 1 and image i; open symbols). For the method of ORC calculation, see text.

Discussion

The calculation of biodiffusion coefficient (D b ) using the ''classical'' method allowed us to confirm the very high impact on sediment reworking of the dominant dense A. filiformis population (Table 3) from a specific station of southeastern Deep Trench of Skagerrak (S5).

The images provided by the optical system have made possible visualization of particle displacements inside the sedimentary column due to the activities of this population. It was anticipated based on results from the classical luminophore method (see above) that bioturbation by this species was of the biodiffusor functional type [START_REF] Franc ßois | A new model of bioturbation for a functional approach to sediment reworking resulting from macrobenthic communities[END_REF]. Instead, it appeared that particles were displaced in a nonlocal mixing pattern sensu [START_REF] Boudreau | Mathematics of tracer mixing in sediments: II. Nonlocal mixing and biological conveyorbelt phenomena[END_REF], i.e. material is transported between points that are not adjacent to each other. Nevertheless, a general trend leading to a homogenization inside the sedimentary column with time was suspected. Indeed, this homogenization pattern is a common trend with time whatever the particular mode of sediment mixing, as illustrated by other specific models (not based on Berner's equation; e.g. [START_REF] Boudreau | Diffusion in a lattice-automation model of bioturbation by small deposit feeders[END_REF]. Those localized movements could be easily linked to the construction of disc chambers and maintenance of burrows by the organisms, as very recently visualized and quantified by [START_REF] Solan | Observation and quantification of in situ animal -sediment relations using timelapse sediment profile imagery (t-SPI)[END_REF].

Our goal was not only to get some images of the 2-D repartition of luminophores in the sediments but also to use this information to quantify sediment reworking by infaunal activities. By analogy with the ''classical'' luminophore method, the calculated amount of tracer was projected in horizontal layers to get a 1-D repartition. Using the obtained profiles that were of diffusive-like shape, biodiffusion coefficients (D b ) were calculated using the Cochran model (1985). However, extremely high D b were then obtained, e.g. 274 cm À 2 year À 1 (data not shown). This was explained by (i) the reduced information contained in the 2-D compared to the 3-D sampling (classical method), and (ii) initial conditions were not appropriate for the Cochran model, i.e. at the start of experiment, t = 0, all the tracers have to be homogeneously located at the sediment surface (C(0;0) = C 0 ; C(z;0) = 0; z: depth).

Obviously, while looking at our images, these requested initial conditions were not fulfilled even when using a 1-cm-thick layer virtual slicing for calculation (e.g. as shown in Fig. 3). Then, in order to use our precise two-dimensional data, we developed a new particle reworking quantification (ORC, optical reworking coefficient; DQ cm À 2 h À 1 ) expressing the total variation of luminophores per surface and time units.

When taking into account the four faces of the rectangular sampling cores, the new 2-D optical method and calculation of the ORC allowed us to efficiently compare the sediment reworking occurring in different samples from the same location. Moreover, no specific initial conditions have to be satisfied in order to calculate the ORC, the coefficient being independent of both the tracer repartition (with potential limitation limit, cf. below) and the size of the quantification cells used (as a slope of a linear relationships). Furthermore, the ORC coefficient appears to be very flexible in use being independent of the time interval between the images. This, however, may have some limitation due to the actual experimental setup that we used (small surface cores), which may have induced a change in the organism behavior inside the sediments after a certain period of incubation. Also, the homogenization of the luminophores repartition with time could reduce the visualization of the reworking process while increasing the quantity of undetected events, i.e. the more homogeneous is the luminophore repartition, the more frequent is the occurrence of an exchange between two fluorescent tracers.

Since the 2-D approach cannot provide biodiffusion coefficients, it cannot replace the 1-D approach, but should rather be regarded as a supplementary technique to study the sediment bio-reworking. The classical method, which is destructive, quantifies (D b ) the global movements of particles inside the sedimentary column but cannot be used for small time and spatial scale measurements. On the other hand, the new optical nondestructive technique allows the quantification of the 2-D particle displacements (ORC) based on high spatial and temporal resolutions. Such high resolution is required to precisely study infauna behavioral ecology as it has been done here with an A. filiformis population. Moreover, contrary to the biodiffusion coefficient, the ORC is not related to a specific type of particle movements. Physically, it corresponds to a quantity of particle displaced per surface unit and per time unit. In the frame of bioturbation studies, it can be then useful to compare sediment reworking between locations whatever the bioturbating organisms (functional groups) or communities involved, whereas D b can be only estimated from diffusive-like particle transport.

The next step in the development of this new method will be to test it on several other types of sediments with other monospecific faunal communities (functional groups) and then to more complex multispecific communities. Based on such studies, the linkage between the 1-D (D b ) and the 2-D (ORC) methods should be, if it exists, established. Moreover, this 2-D optical quantification of sediment reworking is very promising while able to provide a important basis information about the bioturbation process when combined with the recently developments of solute 2-D distribution quantification.

Fig. 1 .

 1 Fig. 1. Map showing the stations visited during the August 1999 and May 2001 (only S5) cruises.

Fig. 3 .

 3 Fig. 3. A. filiformis (S5) and two-dimensional optical luminophore method. Sequence of images (i.e. A to F) showing the 2-D distribution of luminophores (visualized in white color) in the sediments with time. The time increment of each picture, since picture A, is indicated below the images. Scale bars = 1 cm. To improve image visualization, binary but not grayscale images are presented.

Table 1

 1 Sediment reworking parameters, biodiffusion coefficient (D b ) and burying depth, measured during the 1999 cruise with the ''classical luminophore method''

	Station

Table 3

 3 Station S5 abundance per m 2 of the dominant species in the sediments sampled in August 1999 (mean F S.D.,

	n = 3)
	Species

Acknowledgements

We would like to thank Kristina Sundba ¨ck (Marine Botany, Go ¨teborg University, Sweden) who hosted us to perform the luminophore counts on her stereo UV microscope, and Eric Gilbert (BCEOM, Montpellier, France) for his advice on image analysis. Thanks are also due to the crews of N/V Arne Tiselius for their help on the field, to the Benthic Dynamics Conference 2002 participants, especially Yoko Furukawa, Joe Germano, Robert Kennedy and Martin Solan, for their positive feedback, and to two anonymous reviewers for their constructive comments. This work was supported by a mise a `disposition fellowships from the Centre National de la Recherche Scientifique to Franck Gilbert. Financial support to Stefan Hulth was obtained from the Swedish Natural Science Research Council (VR), the Swedish Engineering Foundation (CF), the Technology Link Foundation in Go ¨teborg (TBSG) and the Foundation for Strategic Environmental Research (MISTRA). Nereis Park contribution no. 003. [RW]


Keywords:

This paper demonstrates that Magnitude Squared Incoherence (MSI) analysis is efficient to localize hot spots, i.e. points at which focused Electromagnetic (EM) Analyses can be applied with success. It is also demonstrated that MSI may be applied to enhance Differential EM Analyses (DEMA) based on Difference of Means (DoM).

I. INTRODUCTION

A mong Side-channel attacks some exploit the timing behavior of I [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF], while others exploit the power consumption [START_REF] Kocher | Differential Power Analysis[END_REF] or the EM emissions [START_REF] Gandolfi | Electromagnetic Analysis: Concrete Results[END_REF][START_REF] Peeters | Power and electromagnetic analysis: Improved model consequences and comparisons, Integration[END_REF][START_REF] Agrawal | The EM Side-Channel(s)[END_REF].

EM side channel efficiency is due to the inner properties of EM emissions. Their ability to propagate through different materials is the most striking one. Indeed, it allows attackers targeting the bounded hardware area integrating the cryptomodule under attack or part of it. This is all the more interesting because it allows getting round global hardware countermeasures against power analysis [START_REF] Shamir | Protecting Smart Cards from Passive Power Analysis with Detached Power Supplies[END_REF] by focusing the analysis on reduced silicon die areas thanks to tiny probes. Moreover, replacing the large probe of Fig. 1a by the tiny one shown Fig. 1b allows dividing the number of EM traces required to disclose the key of standard iterative DES (Data Encryption Standard) mapped into a FPGA by 100 to 200. However, this requires localizing beforehand the crypto-module that may occupy only a small fraction of the device area.

However, focusing EM attacks, using small sensors, requires localizing leaking spots to overcome the quadratic increase (with the square of package / probe sizes) of the number of points to be attacked using either Differential EM Analyses (DEMA) [START_REF] Gandolfi | Electromagnetic Analysis: Concrete Results[END_REF][START_REF] Peeters | Power and electromagnetic analysis: Improved model consequences and comparisons, Integration[END_REF] or Correlation EM Analyses (CEMA). Within this context, the contribution of this paper is twofold. First, a localization technique based on Spectral Coherence analyses is introduced. It allows finding positions where EM attacks might be successful with a reduced set of traces. The technique, called Weighted Global Magnitude Squared Incoherence (WGMSI) analysis, has several interesting properties. Firstly, it requires only few EM measurements to Manuscript received September 14, 2010. All Authors are with the Laboratory of Informatics Robotics and Microelectronics of Montpellier, France (email: firstname.name@lirmm.fr) be efficiently applied. Thirdly, this non invasive and contactless technique can be applied even if the circuit under attack is encapsulated.

The second contribution is a Spectral Coherence based technique allowing performing successful DEMA with a reduced set of traces. This technique, which accelerates the convergence of DEMA based on DoM test, is denoted by Differential Global Magnitude Square Incoherence (DGMSI) analysis afterward.

The paper is organized as follows. Section II highlights the localization problem and gives details about the proposed localization technique. This section ends with experimental results demonstrating that WGMSI is efficient to localize the less robust positions against EM attacks. Section III highlights how Spectral Coherence can be further used to enhance DEMA based on DoM by reducing the number of traces to be collected and processed to disclose the secret key. Note that DoM is not the sole statistical test that can be used during an EM attack. Indeed there are many other tests; a classification of test is given in [START_REF] Standaert | Partition vs. Comparison Side-Channel Distinguishers: an Empirical Evaluation of Statistical Tests for Univariate Side-Channel Attacks against Two CMOS Devices[END_REF]. Conclusion is drawn in section IV.

II. WEIGHTED GLOBAL MAGNITUDE SQUARED INCOHERENCE

Performing a DEMA, requires collecting a large number of traces. It is thus time consuming even if the analysis is done at a single position above the device, with a large probe.

The situation becomes unpractical when tiny sensors are used by an attacker or an evaluator aiming at demonstrating the robustness of a design against EM analyses.

A. Problem

Considering an iterative DES mapped onto a FPGA, one may plan collecting, at 784 (28×28) positions above and around the Spartan3 core (displacement step of 250µm), 20 000 EM traces, using a 500µm diameter magnetic loop, to determine if EM emissions may be exploited or not. However, this would result, in collecting traces during 24 days at first order and with the setup presented in annex A.

Enhancing Electromagnetic Analysis Using

Magnitude Squared Incoherence This leads to the question: how efficiently and quickly position tiny probes? One may think, to solve this problem, by performing a EM Near Field Scan (NFS) to localize the cryptographic block. However this is often insufficient.

As an illustration, Fig. 2c shows a peak to peak cartography (obtained with the setup introduced in annex A) of the EM emissions measured at several coordinates above the FPGA package during a DES ciphering. Fig. 2b is an X-Ray photography of the package containing the circuit under attack. Finally, Fig. 2a discloses the floorplan of the circuit, running at 50MHz and integrating a DES module, a finite state machine and a RS232 interface for communication purpose.

As shown, it appears impossible to correlate Fig. 2a with Fig. 2c even if the die area is known thanks to the X-Ray photography. It appears all the more difficult to localize the DES module and thus decide where to position the magnetic sensor to perform a successful focused DEMA or CEMA to avoid potential global hardware countermeasures.

B. Basics of Global Magnitude Squared Incoherence

DEMA exploits the data dependent behavior of the EM emissions radiated by circuits during cryptographic operations. EM emissions are generated by flows of electrical charges through the wires connecting logic gates but also trough wires supplying the circuit [START_REF] Ordas | Near-Field Mapping System to Scan in Time Domain the Magnetic Emissions of Integrated Circuits[END_REF]. Since the switching of gates generates a current flow through the circuit interconnect, we may conclude that gates generate some data dependent EM emissions at different points above the circuit according to the power distribution network. These data dependent behaviors is exploited by statistical means [START_REF] Gandolfi | Electromagnetic Analysis: Concrete Results[END_REF][START_REF] Peeters | Power and electromagnetic analysis: Improved model consequences and comparisons, Integration[END_REF], to retrieve the secret key.

Even if the timing and power characteristics of CMOS gates are known, it is difficult to deduce any characteristic about the EM emissions generated of actual IC due to the complexity of their power distribution grid. Thus, the only conclusion we may consider is that gates generate some EM perturbations i.e. generate some data dependent harmonics located somewhere in the whole EM emission spectrum.

Within this context, the proposed technique allows disclosing the data dependent behavior of EM emissions in the frequency domain without making any assumption on the EM emission characteristics. It is based on spectral incoherence analysis of two time domain signals as detailed below. The only observation on which is based the method is: considering two successive hardware operations, we are sure that some gates switch during one computation and do not switch during the other, while some gates switch during both operations. This leads to the following intuitive conclusion that guides the development of our proposal: between two cryptographic operations some characteristics of the EM emissions remain constant (coherent) from one operation to another, while some characteristics completely change (are incoherent). Such a data dependent behavior is disclosed by the WGMSI technique.

The Magnitude Squared Coherence (MSC) between two signals w1(t) and w2(t) is a real-valued function of frequency with values between 0 and 1. It is defined by:

) ( ) ( ) ( ) ( 2 , 2 1 , 1 2 2 , 1 2 , 1 f P f P f P f MSC w w w w w w w w ⋅ = (1) ) ( 1 ) ( 2 , 1 2 , 1 f MSC f MSI w w w w - = (2)
where P w1w1 (f), P w2w2 (f) are the power spectral density of w1(t), w2(t), and P w1w2 (f) is the cross power spectral density of w1(t) and w2(t). At a given frequency f, a MSC(f) value of 1 indicates that the two spectra are exactly the same while, a value of 0 means that the spectra are different i.e. incoherent. Alternatively, one may compute the Magnitude Squared Incoherence MSI(f) [START_REF] Kocher | Differential Power Analysis[END_REF]. This criterion has also its values between 0 and 1 but indicates rigorously the contrary of (1).

Considering the whole spectra of two time domain signals, one may compute the WGMSI coefficients between them according to (3) that consider the signal w2(t) as a reference.

∑ ∈ ∈ ⋅ = BW f w BW f w w w f A Max f A f MSI WGMSI )) ( ( ) ( ) ( 2 2 2 , 1 (3)
where nf is the number of frequency values at which the MSI(f) coefficients are computed, BW is the considered frequency bandwidth and A w2 (f) is the power spectrum amplitude at the frequency f. WGMSI has values between 0 and 1. A high value indicates that w1(t) and w2(t) have perfectly incoherent spectra, while a low value indicates the contrary. Note that the second term of (3) is a key term. Indeed, it weights MSI(f) values such that incoherent and high amplitude harmonics have more impact on the final WGMSI value than incoherent but low amplitude harmonics. This reduces significantly the impact of noise.

To illustrate these definitions, 5 time domain EM traces were acquired during 5 different data processing of a DES. These traces (Fig. 3) have been collected with a 500µm diameter probe placed respectively above a DES (Fig. 3a) and above a clock wire (Fig. 3b). As a result, one may expect that curves Fig. 3a are data dependent traces while, waveforms Fig. 3b are completely data independent. To validate this assumption, MSC(f) were computed. Fig. 3c gives the MSC(f) evolution with respect to frequency for both data dependent and fully data independent traces. As shown, the MSC(f) values obtained considering traces collected above some clock nets have, as expected, values closer to 1 over a wider frequency range than the MSC(f) values computed with traces collected above the DES, validating the above discussion. The obtained MSC(f) values were gathered to compute the WGMSI coefficients. As expected, WGMSI values (Table 1) corresponding to acquisitions above clock nets are two magnitude order lower than those acquired above the DES. Considering these results, one may assume that the WGMSI criterion appears efficient to differentiate a data dependent behavior from a data independent one and may be used during a magnetic NFS to localize area with data dependent EM emissions expected to be hot spots for EM analysis.

C. Coupling WGMSI and EM Near Field Scan

Coupling WGMSI with NFS system to localize points characterized by data-dependent EM emissions is straightforward. The idea is to collect at each (X,Y) coordinates above the IC at least two different time domain traces of the EM field corresponding to two different data processing. Then, WGMSI values are computed for all (X,Y) positions to obtain a WGMSI cartography revealing positions with data-dependent EM emissions. Note that computing WGMSI values for more than 2 data and averaging the results is not mandatory but leads to better results in practice.

D. WGMSI vs EM cartographies

To validate the effectiveness of the WGMSI analysis, two kinds of validation were performed. One aimed at correlating the obtained WGMSI cartographies with design floorplans while the 2 nd aimed at demonstrating that spots with the highest WGMSI values are good candidates for DEMA and CEMA. Results related to the 1 st validation step being available in [START_REF] Dehbaoui | Enhancing Electromagnetic Attacks using Spectral Coherence based Cartography[END_REF] we give here only the results related to the 2 nd validation step. Note however, that the 1 st validation step allowed us concluding that WGMSI technique is an efficient technique to localize a cryptographic system, and its Power/Ground (P/G) network, within noise.

Note finally that both validation steps have been carried out with the experimental setup described in annex A and were performed considering a design mapped into a FPGA circuit and more precisely a Spartan3 board. The mapped design (50MHz) integrated: a RS232 to communicate with the PC, a finite state machine that manages the communications and the behavior of the chip, and a DES.

To evaluate the efficiency of WGMSI cartographies in disclosing hot spots for EMA, we scanned a 7mm × 7mm area of the package centered on the Spartan core with a 1mm probe displacement step. The core size was estimated to be roughly of 5×5mm² thanks to an X-ray photography, while the package has a size of 17× 17mm².

At each of the 64 resulting positions, 5000 averaged traces (20 trials for each PTI) were collected corresponding to the ciphering of 5000 different PTI. Note that the same 5000 PTI were processed at each position in the same order. Data acquisition achieved, WGMSI cartographies were computed considering two frequency bandwidths, [40-, 950MHz] and , and EM traces corresponding to the whole DES course. The 1 st bandwidth is the full bandwidth of our acquisition chain, while the 2 nd one has been chosen to keep only the emissions of the P/G network according to [START_REF] Ordas | Near-Field Mapping System to Scan in Time Domain the Magnetic Emissions of Integrated Circuits[END_REF].

Tables II andIII give the WGMSI values in both cases after normalization done so that the maximum WGMSI value corresponds to 100% and the minimum one to 0%. Note that black cells point out positions with WGMSI values lower than 20%. As shown, the DES and part of P/G rails in which flows its switching current generate data dependent EM emissions mainly in left side of the cartographies.

In a 2 nd stage, 64 DEMA and CEMA, considering the Hamming Distance model were performed; both attacks targeted the last round of the DES. To compare the results obtained at different positions two criteria were considered. Note we performed CEMA using the Pearson's correlation. Its evaluation resumes in detecting and storing, while the number of processed traces increases, the beginning of the 1 st sequence of 100 successive right guesses of the key. Comparing these tables highlights of several results. First, Tables II and III indicate that the DES and part of the overall P/G network supplying it are localized mainly on the left of the core as expected from the floorplan.

Second, Tables II to VII show that WGMSI and MTDwS cartographies have similarities. Indeed, most positions with high MTDwS values are mainly on the right part of the cartographies as positions with low WGMSI values. To better support this observation, correlations between WGMSI cartographies and, all other ones have been computed. Table VIII gives the results. Note that to compute the correlations involving MTDwS, we considered MTDwS equal to 100 for positions where attacks were unsuccessful. As shown, there are correlations between these cartographies. This demonstrates the interest of WGMSI to localize hot spots. Note that higher correlation values are obtained for the frequency bandwidth [40MHz -200MHz]. According to [START_REF] Ordas | Near-Field Mapping System to Scan in Time Domain the Magnetic Emissions of Integrated Circuits[END_REF], this means that the most leaking points are part of the P/G network.

] [ ] ( ) [ ] ( ) ∑ ∑ ∑ ∑ = = = = - ⋅ - - ⋅ = N i s i N i i s i N i s i N i i s i Ks ) ,K D(PTI j T ) ,K D(PTI ) ,K D(PTI j T ) ,K D(PTI j ∆ 1 1 1 1 1 1 (4)
In ( 4), extracted from [START_REF] Bevan | Ways to Enhance Differential Power Analysis[END_REF], ∆ Ks [j] is the j th sample of the DoM, N is the number of EM traces used, PTI i is the i th plaintext, T i [j] is the j th sample of associated EM trace and D the selection function returning the value of the targeted bit according to the PTI and the sub-key guess.

According to [START_REF] Kocher | Timing Attacks on Implementations of Diffie-Hellman, RSA, DSS, and Other Systems[END_REF], if the sub-key guess is correct the right and left hand terms of (4) correspond to the averaged traces characterized by an effective targeted bit value of 1 and 0, respectively. Thus, the difference of these terms ∆ Ks [j] shows a bounce as illustrated by Fig. 4b.

Contrarily, if the sub-key guess is wrong, the left and right terms are undistinguishable, and a smaller bounce appears. Indeed, the selection function D fails in sorting traces according to the targeted bit and these terms are expected to be less incoherent than for the right guess.

Considering this, one may compute, for each sub-key guess Ks, the Global Magnitude Squared Incoherence, GMSI Ks , between these two means curves. Then, rather than computing the DoM, one may compute the DoM weighted by the GMSI Ks [START_REF] Agrawal | The EM Side-Channel(s)[END_REF], and search the sample with the highest amplitude:

[ ] [ ] ∑ ∈ = ⋅ ∆ = ∆ BW f Ks Ks Ks Ks nf f MSI GMSI GMSI j j ) ( * (5)
Let us denote by DGMSI, a DEMA attack considering the four output bits of Sboxes [START_REF] Bevan | Ways to Enhance Differential Power Analysis[END_REF] and performed using the proposed weighting strategy.

To evaluate its efficiency, we computed MTDwS and PRG cartographies using only 1000 PTI rather than 5000. Tables IX and X give the results. As shown, a DGMSI performed with 1000 PTI gave results quite similar to that of a DEMA performed with 5000 PTI. One may conclude that DGMSI allows reducing by ∼5, wrt DEMA, the number of PTIs required to disclose the secret key. This demonstrates the interest of MSI to enhance EM analyses. One may wonder why DGMSI allows reducing MTDwS and increasing the PRG. Practically, the obtained enhancement is explained by the reduction of the set of key candidates during the key search by removing keys related to coherent means in the DoM. Fig. 4 shows the DoM obtained with DGMSI and DEMA after the processing of the 1000 PTI. As shown, the right key appears more clearly with DGMSI than for DEMA: the margin between the right and wrong keys is 2 times greater with DGMSI than with DEMA. This is due to the low values of GMSI Ks that flatten the DoM of wrong keys. In the previous section, the weighting strategy of DoM has been demonstrated efficient to reduce the number of EM traces required to disclose the secret key. However, during this experiment, we considered averaged EM traces. One may wonder if results hold in presence of noise. 5000 one shot EM traces, were thus acquired and processed using DEMA, CEMA and DGMSI. Note that 10 successive attacks were done; each one differing from the other by the order considered to process data. Table XI gives the results. In this table, µ and σ denotes the mean and the standard deviation obtained considering the 10 processing orders.

As shown, with one shot EM traces, DGMSI allows disclosing the right key with less PTI than with DEMA. Indeed, only 43% of the traces are processed in average to reach for the first time the sequence of 100 successive right guesses of the full key with DGMSI while 68% and 73% are needed with DEMA and CEMA respectively. Moreover the PRG of DGMSI is nearly twice the ones of DEMA and CEMA.

IV. CONCLUSION

The interest of Magnitude Squared Incoherence processing technique for EM analyses has been demonstrated through two applications. One tackles the identification of hot spots to be attacked using tiny probes while the other is an improvement of DEMA based on DoM.

Annex A All results were obtained with a measurement platform composed of:

-an oscilloscope, with a 3.5 GHz bandwidth, to sample at 40 GS/s the time domain evolutions of the EM signals, -a low noise 63 dB amplifier with a 1GHz bandwidth, -a magnetic loop with a 500µm diameter, and a bandwidth greater than 1GHz, -a motorized stage allowing positioning along X, Y and Z axes the probe with a resolution of 10µm, -a PC to control the whole measurement setup, i.e. to provide data to the DES (50 MHz) through an RS232 module and store the EM traces collected by the scope.

AFig. 1 .

 1 Fig. 1. (a) global EM analysis and (b) focused EM analysis

  die

Fig. 2 .

 2 Fig. 2. (a) Design floorplan (b) X-Ray photography of the package (c) measured Peak to Peak EM cartography above the package.

  The 1 st criterion was the Measurements To Disclosure with Stability (MTDwS) adopted during the DPA contest 2009 [8].

Fig. 3 .

 3 Fig. 3. EM traces collected (a) above the DES (b) above the clock network and (c) MSC(f) versus frequency

[

Fig. 4 .

 4 Fig. 4. DoM obtained with (a) DGMSI and (b) DEMA

TABLE I .

 I WGMSI VALUES

TABLE II .

 II NORMALIZED (%) WGMSI CARTOGRAPHY / [40MHZ, 950MHZ]

	Y/X	1	2	3	4	5	6	7	8
	8	33	29	14	5	3	1	1	1
	7	53	35	24	10	4	1	1	1
	6	85	57	50	33	13	1	1	6
	5	90	62	64	63	52	13	20	11
	4	83	48	30	16	13	1	1	2
	3	100	37	20	12	8	2	1	1
	2	75	24	17	12	6	3	2	2
	1	40	24	20	12	7	3	2	3
		TABLE III. NORMALIZED (%) WGMSI / [40MHZ, 200MHZ]
	Y/X	1	2	3	4	5	6	7	8
	8	9	12	14	9	4	1	1	1
	7	8	12	23	22	9	1	1	2
	6	11	11	26	27	27	1	2	14
	5	30	12	33	38	41	6	100	32
	4	39	10	21	20	22	2	4	8
	3	57	30	8	11	4	4	1	4
	2	36	21	5	6	8	6	2	1
	1	27	5	1	2	2	3	5	3

  Table IV and V show normalized MTDwS values for DEMA and CEMA. The normalization has been done so that a value of 60% means that 60% of EM traces have been processed to reach the 1 st sequence of 100 successive right guesses. Black cells indicate positions at which the MTDwS is greater than 50% or positions at which the attack failed.Finally, the 2 nd criterion is the Percentage of Right Guesses (PRG) obtained after the processing of 5000 traces. Tables VI and VII show results. Black cells point out positions with a PRG value lower than 20%.

TABLE IV .

 IV NORMALIZED MTDWS (%) VALUES FOR DEMA (5000 PTIS)

	Y/X	1	2	3	4	5	6	7	8
	8	84	58	36	62	fail	fail	fail	fail
	7	62	34	98	86	88	fail	fail	fail
	6	98	94	86	48	48	fail	fail	52
	5	58	82	80	58	34	48	40	34
	4	36	56	70	fail	48	fail	fail	80
	3	96	94	54	52	98	fail	fail	fail
	2	50	74	fail	48	50	fail	84	fail
	1	80	fail	fail	70	42	fail	fail	fail
	TABLE V. NORMALIZED MTDWS (%) VALUES FOR CEMA (5000 PTIS)
	Y/X	1	2	3	4	5	6	7	8
	8	96	36	38	40	60	fail	fail	72
	7	50	35	35	32	39	61	41	81
	6	31	38	19	20	38	95	23	35
	5	36	53	21	9	38	29	25	12
	4	28	79	13	15	41	35	34	21
	3	7	40	18	18	16	18	fail	32
	2	49	43	63	17	22	24	33	fail
	1	23	71	fail	69	30	32	39	fail
	TABLE VI. NORMALIZED PRG (%) VALUES FOR DEMA (5000 PTIS)
	Y/X	1	2	3	4	5	6	7	8
	8	40	60	65	39	0	0	0	0
	7	38	71	33	40	54	0	0	0
	6	33	57	10	69	57	0	0	49
	5	11	36	17	35	66	53	64	66
	4	66	44	13	0	61	0	0	42
	3	65	17	53	60	27	0	0	1
	2	51	6	0	69	51	0	6	0
	1	46	0	0	19	67	0	2	0
	TABLE VII. NORMALIZED PRG (%) VALUES FOR CEMA (5000 PTIS)
	Y/X	1	2	3	4	5	6	7	8
	8	40	71	64	61	50	5	0	43
	7	54	70	75	75	72	39	64	32
	6	77	72	88	85	62	52	78	65
	5	65	62	85	92	63	75	80	89
	4	77	41	91	90	61	69	67	85
	3	80	62	85	86	84	90	0	82
	2	55	67	49	84	85	87	69	0
	1	78	21	0	40	70	69	64	0

TABLE VIII .

 VIII CORRELATIONS BETWEEN WGMSI AND MTDWS, PRG VALUES

		Correlation	MTDwS	PRG
	DEMA	WGMSI 40 -950 MHz	-0.33	0.41
		WGMSI 40 -200 MHz	-0.46	0.49
	CEMA	WGMSI 40 -950 MHz	-0.20	0.20
		WGMSI 40 -200 MHz	-0.40	0.38
	III. DIFFERENTIAL GLOBAL MAGNITUDE SQUARED
		INCOHERENCE
	If the above results have demonstrated the interest of the
	WGMSI, they also suggest using MSI to enhance DEMA.
	Attacking a DES sub-key by DEMA consists in computing
	64 DoM according to (4) and in identifying the one having the
	sample with the highest amplitude.

TABLE IX .

 IX NORMALIZED MTDWS (%) VALUES FOR DGMSI (1000 PTIS)

	Y/X	1	2	3	4	5	6	7	8
	8	91	96	48	92	fail	fail	fail	fail
	7	fail	71	22	24	60	fail	fail	fail
	6	fail	24	29	24	40	fail	fail	25
	5	28	60	28	28	28	41	21	26
	4	25	fail	28	29	24	fail	fail	60
	3	21	43	65	36	26	fail	fail	fail
	2	27	39	69	27	42	65	fail	fail
	1	41	58	fail	fail	63	fail	fail	fail
	TABLE X. NORMALIZED PRG (%) VALUES FOR DGMSI (1000 PTIS)
	Y/X	1	2	3	4	5	6	7	8
	8	11	12	63	30	0	0	0	0
	7	10	72	80	78	46	0	0	0
	6	16	76	73	77	69	0	16	75
	5	74	62	72	78	76	59	79	75
	4	76	9	75	71	76	0	0	33
	3	79	60	35	76	74	18	0	0
	2	73	68	38	86	62	37	0	0
	1	59	44	0	1	42	10	0	0

TABLE XI .

 XI COMPARING DGMSI WITH DEMA AND CEMA (5000 PTI)

		DGMSI		DEMA	CEMA
		MTDwS	PRG	MTDwS	PRG	MTDwS	PRG
	µ	43%	58%	68%	33%	73%	26%
	σ σ σ σ	11%	10%	17%	17%	19%	20%


F Gilbert
email: gilbert@com.univ-mrs.fr

G Desrosiers

S Hulth

G Stora

Comparison between the Nereis diversicolor and Nereis virens marine worms in the transformation of ingested hydrocarbons

A feeding experiment was conducted on the marine worm Nereis diversicolor to compare the fate of a hydrocarbon mixture during the gut passage in this species with the hydrocarbon breakdown process demonstrated for Nereis virens. Hydrocarbon dissolution/solubilization processes in the gut of N. diversicolor were found to have similar qualitative and quantitative importance in the hydrocarbon transformation as those observed in N. virens.

Macrobenthic worm species are known to play an important role in the vertical redistribution of hydrocarbons in the sedimentary column and subsequent degradation by microbes [START_REF] Lee | Biological processes a¡ecting the distribution of pollutants in marine sediments. Part II. Biodeposition and bioturbation[END_REF][START_REF] Mcelroy | In£uence of mode of exposure and the presence of a tubiculous polychaete on the fate of Benz(a)anthracene in the benthos[END_REF]. In areas of intensive reworking, ingestion of hydrocarbons adsorbed to sediment particles may induce direct hydrocarbon transformation during the transit. [START_REF] Mayer | Bioavailability of sedimentary contaminants subject to deposit-feeder digestion[END_REF] discovered a potentially powerful dissolving process for organic matter related to the production of surfactants inside the gut of the deposit-feeders. Most recently, this process was demonstrated and quanti¢ed for Nereis virens (Sars) species [START_REF] Bock | Digestive plasticity of the marine benthic omnivore Nereis virens[END_REF][START_REF] Gilbert | Alteration and release of aliphatic compounds by the polychaete Nereis virens (Sars) experimentally fed with hydrocarbons[END_REF]. The aim of our present work was to investigate if Nereis diversicolor (O.F. Mu« ller) has a similar capacity to transform ingested hydrocarbons as Nereis virens and thus plays an important role in the transport and transformation of hydrocarbons in bioturbated surface deposits.

Eighteen adult N. diversicolor were collected in the Carteau Cove (Gulf of Fos, Mediterranean Sea, France) in February 2000. A similar experimental design was employed as that previously used for N. virens [START_REF] Gilbert | Alteration and release of aliphatic compounds by the polychaete Nereis virens (Sars) experimentally fed with hydrocarbons[END_REF] adapted to N. diversicolor (animal size, period of faeces evacuation). After a 7-d starvation period, each worm was introduced in a horizontal glass tube. Tubes were placed in experimental containers (three worms per container) ¢lled with ¢ltered and aerated seawater under constant temperature (158C) and L:D cycle (12 h:12 h). Once a day, the polychaetes were individually fed with a portion of fresh algae Enteromorpha spp. (9.4 AE1.1mg; mean AESD; N216) anaesthetized with a mixture of alkanes (n-C 14 , n-C 15 , n-C 17 , n-C 18 , n-C 22 , n-C 23 , n-C 28 and n-C 29 , pristane and androstane: 2.4 AE0.3 mg; mean AESD; N216). After four, eight, 12 and 24 h, the faeces were sampled in each container, weighed, and frozen. The sum of faeces collected during the 24 h period after worm feeding was used for the analysis of the aliphatic hydrocarbons. Eventual food residues were removed, and the sea water replaced prior to feeding. The worms were fed over a period of 12 d. After ten, 11 and 12 d, the worms in one container were sacri¢ced, and their gut tissues isolated, weighed and frozen before analysis of aliphatic hydrocarbon content. In an additional container without worms, three food rations were separately subjected to the same experimental conditions for 24 h and served as controls. The aliphatic hydrocarbons were extracted from the homogenized algae, worm faeces and gut tissues and analysed by GC/ MS [START_REF] Gilbert | Alteration and release of aliphatic compounds by the polychaete Nereis virens (Sars) experimentally fed with hydrocarbons[END_REF]. Di¡erences between species were studied using a one-way ANOVA after transformation of the heteroscedastic data.

Hydrocarbon content in the faeces of N. diversicolor were found to be more or less constant throughout the duration of the experiment (Table 1). A similar observation was made in N. virens [START_REF] Gilbert | Alteration and release of aliphatic compounds by the polychaete Nereis virens (Sars) experimentally fed with hydrocarbons[END_REF]. In order to give a readable presentation of the data, three randomly chosen sets of data were used as triplicates (corresponding to days eight, nine and 11) to compare results between the two species.

The quantitative distributions of hydrocarbon tracer in gut tissues and faeces were di¡erent for the two species (Table 1, Figure 1). After ingestion of the algal material, the gut tissue content were not signi¢cantly di¡erent between the two species (P0.24, Table 1). However, the percentage of tracer in faeces with regard to the initial food portion content was signi¢cantly higher (P40.05) for N. diversicolor than for N. virens. Only 36.8 and 51.1% of the ingested hydrocarbon tracer was found associated with worm gut tissues and solid material (faeces), for N. virens and N. diversicolor, respectively, suggesting greater rejection of the tracer by N. virens as non-solid material.

All the compounds of the initial tracer mixture were present in the control, N. diversicolor and N. virens faeces (Figure 1). In both the control and the faeces there was a reduction of the alkane content, which was most pronounced for the low molecular weight n-alkane (LALK; from n-C 14 to n-C 18 ) than the heavy molecular weight n __ -alkanes (HALK; from n-C 22 to n-C 29 ). An increased loss of the low molecular weight compounds was apparent in the faeces, whereas the di¡erences measured for the highest molecular weight n-alkanes had decreased. This resulted in a global decrease of the LALK/HALK ratios in the faeces compared to the control (Figure 1). Moreover, this tendency (P0.26) seems to be more important for N. virens than for N. diversicolor.

As previously demonstrated in N. virens [START_REF] Gilbert | Alteration and release of aliphatic compounds by the polychaete Nereis virens (Sars) experimentally fed with hydrocarbons[END_REF], the ingestion and excretion of the hydrocarbon mixture by N. diversicolor stimulated the loss of hydrocarbons compared to the controls. This may suggest that similar processes, such as increased organic matter dissolution/solubilization due to surfactant production inside the gut [START_REF] Mayer | Digestive environments of benthic macroinvertebrate guts: enzymes, surfactants and dissolved organic matter[END_REF][START_REF] Bock | Digestive plasticity of the marine benthic omnivore Nereis virens[END_REF], are also important for N. diversicolor. However, it seems that subsequent loss of hydrocarbons in the liquid phase (49.1%) was lower for N. diversicolor than for N. virens (63.7%) suggesting, in addition with the hydrocarbon gut and faeces content measured, a reduced alkane absorption e¤ciency by N. diversicolor [START_REF] Ahrens | The role of digestive surfactants in determining bioavailability of sediment-bound hydrophobic organic contaminants to two deposit-feeding polychaetes[END_REF]). The lower digestive potential for N. diversicolor could be attributed to: (i) the di¡erence in body size between worms in the two studies (3.16 0.29 g, mean AESD, N20, [START_REF] Gilbert | Alteration and release of aliphatic compounds by the polychaete Nereis virens (Sars) experimentally fed with hydrocarbons[END_REF], and 0.57 0.05 g, N18 for N. virens and N. diversicolor) as demonstrated within a Nereis succinea population (Ahrens et al., in press); and (ii) the thermal experimental conditions used in the studies. Indeed, a temperature of 158C may have enhanced general activities for the N. virens species sampled in the Saint-Lawrence estuary, Canada [START_REF] Desrosiers | Cycle de de¨veloppement d'une population intertidale de Nereis virens (Polychaeta Nereidae) de l'estuaire maritime du Saint Laurent[END_REF]. It may on the other hand have induced lower activities for the Mediterranean N. diversicolor species that are normally used to more elevated temperatures.

Although some eventual artefacts associated with the experimental set-up, the primary conclusion of our study is that N. diversicolor species induced similar changes to an ingested hydrocarbon mixture than were previously demonstrated for N. virens. These two worms may thus signi¢cantly in£uence the fate and distribution patterns of hydrocarbons in surface sediments.

Table 1. Percentage of the hydrocarbon mixture recovered in gut tissues and faeces of Nereis virens and Nereis diversicolor compared to the initial hydrocarbon food portion content. Values are mean AESD (N3).

Gut tissues Faeces

Nereis virens 3.1 AE2.6 33.7 AE5.5 Nereis diversicolor 1.2 AE0.7 48.9 AE3.3

Figure 1 .

 1 Figure 1. Percentage (mean AESD; N3) of the various hydrocarbon compounds recovered in the control and the faeces of Nereis virens and Nereis diversicolor (sum of faeces collected during the 24 h period after worm feeding) compared to the initial tracer (food portion) composition. Numbers refer to carbons in n-alkanes (n-C 14 to n-C 29 ); Pr, pristane; Andr, androstane; LALK/HALK, ratio between the low molecular weight n-alkanes (LALK; from n-C 14 to n-C 18 ) and the heavy molecular weight n-alkanes (HALK; from n-C 22 to n-C 29 ) ratio.

We thank an anonymous referee and M.J. Ahrens who helped to improve our manuscript. This work was supported by the GDR 1123 `HYCAR' (France) and the Preem Environment Foundation (Sweden).