Olivier Cappé

Randal Douc

Arnaud Guillin

Jean-Michel Marin

Christian Robert

Adaptive Importance Sampling in General Mixture Classes

Keywords: Importance sampling, Adaptive Monte Carlo, Mixture models, Entropy, Kullback-Leibler divergence, EM algorithm, Population Monte Carlo Échantillonnage préférentiel, mélanges de loi, méthodes de Monte-Carlo adaptatives, entropie

In this paper, we propose an adaptive algorithm that iteratively updates both the weights and component parameters of a mixture importance sampling density so as to optimise the performance of importance sampling, as measured by an entropy criterion. The method, called M-PMC, is shown to be applicable to a wide class of importance sampling densities, which includes in particular mixtures of multivariate Student t distributions. The performance of the proposed scheme is studied on both artificial and real examples, highlighting in particular the benefit of a novel Rao-Blackwellisation device which can be easily incorporated in the updating scheme.

Introduction

In recent years, there has been a renewed interest in using Monte Carlo procedures based on Importance Sampling (abbreviated to IS in the following) for inference tasks. Compared to alternatives such as Markov Chain Monte Carlo methods, the main appeal of IS procedures lies in the possibility of developing parallel implementations, which becomes more and more important with the generalisation of multiple core machines and computer clusters. Importance sampling procedures are also attractive in that they allow for an easy assessment of the Monte Carlo error (provided trustworthy estimates of the variance can be produced). As a consequence, it is therefore easier to construct learning mechanisms in IS settings because of this ability to compare the errors. In many applications, the fact that IS procedures may be tuned-by choosing an appropriate IS density-to minimise the approximation error for a specific function of interest is also crucial. On the other hand, the shortcomings of IS approaches are also well-known, including a poor scaling to highly multidimensional problems and an acute sensitivity to the choice of the IS density combined with the fact that it is impossible to come up with a universally efficient IS density. While there exist a wide variety of solutions in the literature (see, e.g. Robert and Casella, 2004, Chapter 14), this paper concentrates on the construction of adaptive importance sampling schemes in which the IS density is gradually improved based on the outcome of previous Monte Carlo draws.

While the method proposed here can be traced back to authors such as [START_REF] West | Modelling with mixtures[END_REF] or [START_REF] Oh | Integration of multimodal functions by Monte Carlo importance sampling[END_REF], it is closely related to the so-called Population Monte Carlo (henceforth abbreviated to PMC) approach-in the sense of an iterated simulation of importance samples and in opposition to Markov Chain Monte Carlo simulation that only produces a point at a time-introduced by [START_REF] Cappé | Population Monte Carlo[END_REF]. We briefly review the PMC approach, following the exposition of [START_REF] Cappé | Population Monte Carlo[END_REF] and Douc et al. (2007a,b), in order to highlight the differences with the present work. In PMC, a sample (X 1 , . . . , X N ) approximately distributed from π, is repeatedly perturbed stochastically using an arbitrary Markov transition kernel q(x, x ) so as to produce a new sample (X 1 , . . . , X N ). Conducting a resampling step based on the IS weights ω i = π(X i )/q(X i , X i ), it is then possible to produce a new unweighted sample ( X1 , . . . , XN ) that also constitutes an approximation to the target distribution π. Adaptivity in PMC was achieved by considering a transition kernel q consisting of a mixture of fixed transition kernels

q α (x, x ) = D d=1 α d q d (x, x ) , D d=1 α d = 1 , (1)
whose weights α 1 , . . . , α D are tuned adaptively, along the iteration of the PMC algorithm. The adaptation procedure proposed by Douc et al. (2007a), termed D-kernel PMC, aims at minimising the deviance or entropy criterion between the kernel q α and the target π,

E(π, q α ) = E X π [D(π q α (X, •))] , (2)
where D(p q) = log{p(x)/q(x)} p(x)dx denotes the Kullback-Leibler divergence (also called relative entropy), and where the expectation is taken under the target distribution X ∼ π since the kernels q d (x, x ) depend on the starting value x. In the sequel, we refer to the criterion in (2) as the entropy criterion since it is obviously related to the performance measure used in the cross-entropy method of [START_REF] Rubinstein | The Cross-Entropy Method[END_REF]. In [START_REF] Douc | Minimum variance importance sampling via population Monte Carlo[END_REF], a version of this algorithm was developed to minimise the asymptotic variance of the IS procedure, for a specific function of interest, in lieu of the entropy criterion.

A major limitation in the approaches of Douc et al. (2007a,b) is that the proposal kernels q d remain fixed over the iterative process while only the mixture weights α d get improved. In the present contribution, we remove this limitation by extending the framework of Douc et al. (2007a) to allow for the adaption of IS densities of the form

q (α,θ) (x) = D d=1 α d q d (x; θ d ) , (3)
with respect to both the weights α d and the internal parameters θ d of the component densities. The proposed updating mechanism is quite similar to the EM algorithm with the E-step replaced by IS computations. As demonstrated through the example considered in Section 4, this adaptive scheme is applicable to very general families of latent-data IS densities. A possible drawback of adapting the internal parameters θ d of the component densities is that it sometimes raises challenging robustness issues, particularly when (multidimensional) scaling parameters are tuned. We thus propose a Rao-Blackwellisation scheme that empirically appears to be very efficient while inducing a modest additional algorithmic complexity.

Note again that we consider here the generic entropy criterion of Douc et al. (2007a) rather than the function-specific variance minimisation objective of [START_REF] Douc | Minimum variance importance sampling via population Monte Carlo[END_REF]. This choice is motivated by the recognition that in most applications, the IS density is expected to perform well for a range of typical functions of interest rather than for a specific target function h. In addition, the generalisation of the approach of [START_REF] Douc | Minimum variance importance sampling via population Monte Carlo[END_REF] to a class of mixture IS densities that are parameterised by more than the weights remains an open question (see also Section 5). A second remark is that in contrast to the previously cited works and as obvious in equation ( 3), we consider in this paper only "global" independent IS densities. Thus, the proposed scheme is based on genuine iterated importance sampling, contrary to what happens when using more general IS transition kernels as in (1). Obviously, resorting to moves that depend on the current sample is initially attractive because it allows for some local moves as opposed to the global exploration required by independent IS densities. However, the fact that the entropy criterion in (2) is a global measure of fit tends to modify the parameters of each transition kernel depending on its average performance over the whole sample, rather than locally. In addition, structurally imposing a dependence on the points sampled at the previous iteration induces some extra-variability which can be detrimental when more parameters are to be estimated.

The paper is organised as follows: In Section 2, we develop a generic updating scheme for independent IS mixtures (3), establishing that the integrated EM argument of Douc et al. (2007a) remains valid in our setting. Note once again that the integrated EM update mechanism we uncover in this paper is applicable to all missing data representations of the proposal kernel, and not only to finite mixtures. In Section 3, we consider the case of Gaussian mixtures which naturally extend the case of mixtures of Gaussian random walks with fixed covariance structure considered in Douc et al. (2007a,b). In Section 4, we show that the algorithm also applies to mixtures of multivariate t distributions with the continuous scale mixing representation used in [START_REF] Peel | Robust mixture modelling using the t distribution[END_REF]. Section 5 provides some conclusive remarks about the performance of this approach as well as possible extensions.

2 Adapting the Importance Sampling Density

The M-PMC Algorithm

When considering independent mixture IS densities of the form (3), the entropy criterion E defined in (2) reduces to the Kullback-Leibler divergence between the target density π and the mixture q (α,θ) :

E(π, q (α,θ) ) = D(π q (α,θ) ) = log π(x) D d=1 α d q d (x; θ d ) π(x)dx . (4)
As usual in applications of the IS methodology to Bayesian inference, the target density π is known only up to a normalisation constant and we will focus on a self-normalised version of IS that solely requires the availability of an unnormalised version of π [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF]. As a side comment, note that while E(π, q (α,θ) ) is a convex function of the weights α 1 , . . . , α D (Douc et al., 2007a), it generally fails to be so when also optimising with respect to the component parameters θ 1 , . . . , θ D . Given that minimising (4) in (α, θ) is equivalent to maximising log

D d=1 α d q d (x; θ d ) π(x) dx , (5)
we are facing a task that formally resembles standard mixture maximum likelihood estimation but with an integration with respect to π replacing the empirical sum over observations. This analogy suggests that it is possible to maximise the entropy criterion in (4) using an approach based on the principle of the EM algorithm and, in particular, the use of the augmented mixture representation (involving the indicator variables associated with each component of the mixture). Before providing the details of the derivation in Section 2.2, we first state below the proposed adaptive IS algorithm which we will refer to as M-PMC (for Mixture PMC ) in the following. Let (X i,t ) 1≤i≤N and (α t,N , θ t,N ) denote, respectively, the IS sample and the estimated mixture parameters at the t-th iteration of the algorithm.

Algorithm 1. (M-PMC Algorithm) At iteration t, 1. Generate a sample (X i,t ) from the current mixture IS proposal (3) parameterised by (α t,N , θ t,N ) and compute the normalised importance weights

ωi,t = π(X i,t ) D d=1 α t,N d q d (X i,t ; θ t,N d ) N j=1 π(X j,t ) D d=1 α t,N d q d (X j,t ; θ t,N d ) (6)
and the mixture posterior probabilities

ρ d (X i,t ; α t,N , θ t,N ) = α t,N d q d (X i,t ; θ t,N d ) D =1 α t,N q (X i,t ; θ t,N ) , (7)
for i = 1, . . . , N and d = 1, . . . , D.

2. Update the parameters α and θ as

α t+1,N d = N i=1 ωi,t ρ d X i,t ; α t,N , θ t,N , θ t+1,N d = arg max θ d N i=1 ωi,t ρ d X i,t ; α t,N , θ t,N log q d X i,t ; θ t,N d , (8)
for d = 1, . . . , D.
The convergence of the algorithm may be monitored by computing the so-called normalised perplexity exp(H t,N )/N , where H t,N = -N i=1 ωi,t log ωi,t is the Shannon entropy of the normalised IS weights. The normalised perplexity provides an estimate of exp[-E(π, q (α t,N ,θ t,N ) )] and, for sufficiently large N , it is nondecreasing with t.

Detailed Derivation

Integrated Updates

Starting from (5), assume for the moment that integration with respect to π is feasible. In order to update the parameters of the independent IS density (3), we will take advantage of the latent variable structure that underlines the objective function (5). The resulting algorithm-still theoretical at this stage as it involves integration with respect to π-may be interpreted as an integrated EM (Expectation-Maximisation) scheme that we now describe. Let α t = (α t 1 , . . . , α t D ) and θ t = (θ t 1 , . . . , θ t D ) denote, respectively, the mixture weights and the component parameters at the t-th iteration of this integrated EM algorithm.

As usual in mixtures, the latent variable Z is the component indicator, with values in {1, . . . , D} such that the joint density f of x and z satisfies f (z) = α z and f (x|z) = q z (x; θ z ) , which produces (3) as the marginal in x. As in the standard EM algorithm, we can then take advantage of this latent variable representation. Since the joint density of X and Z is α z q z (x; θ x ), the expectation corresponding to the E step of the EM algorithm is the expected complete log-likelihood, namely, at iteration t of our algorithm,

E X π E Z (α t ,θ t ) {log (α Z q Z (X; θ Z )) |X} ,
where the inner expectation is computed under the conditional distribution of Z in the mixture model given the current value (α t , θ t ) of the parameters, i.e.

f (z|x) = α t z q z (x; θ t z ) D d=1 α t d q d (x; θ t d ) ,
while the outer expectation is under the distribution X ∼ π.

The proposed updating mechanism then corresponds to setting the new parameters (α t+1 , θ t+1 ) equal to

(α t+1 , θ t+1 ) = arg max (α,θ) E X π E Z (α t ,θ t ) {log(α Z q Z (X; θ Z ))|X} , (9)
as in the regular EM estimation of the parameters of a mixture, except for the extra expectation over X. It is straightforward to check that the convexity argument used for the EM algorithm also applies in this setup and, hence, that (E(π, q (α t ,θ t ) )) t≥1 is a non-decreasing sequence. Setting

ρ d (X; α, θ) = α d q d (X; θ d ) D =1 α q (X; θ ) ,
the maximisation program in (9) reduces to

α t+1 = arg max α E X π D d=1 ρ d (X; α t , θ t ) log(α d ) , θ t+1 = arg max θ E X π D d=1 ρ d (X; α t , θ t ) log(q d (X; θ d )) ,
where the first maximisation to be carried out under the constraint that

D d=1 α t+1 d = 1. Hence, α t+1 d = E X π ρ d (X; α t , θ t ) , (10)
θ t+1 d = arg max θ d E X π ρ d (X; α t , θ t ) log(q d (X; θ d )) . (11)
As in the regular mixture estimation problem, the resolution of this maximisation program ultimately depends on the shape of the density q d . If q d belongs to an exponential family, it is easy to derive a closed-form solution for (11), which however involves expectations under π. Section 3 provides an illustration of this fact in the Gaussian case, while the non-exponential Student's t case is considered in Section 4.

Approximate Updates

To make the previous algorithm practical, adaptivity must be achieved by updating the parameters based on the previously simulated IS sample. We thus start the algorithm by arbitrarily fixing the mixture parameters (α 1 , θ 1 ) and we then sample from the resulting proposal α 1 d q d (x; θ 1 d ) to obtain our initial sample (X i,1 ) 1≤i≤N , associated with the latent variables (Z i,1 ) 1≤i≤N that indicate from which component of the mixture the corresponding (X i,1 ) 1≤i≤N have been generated. From this stage, we proceed recursively. Starting at iteration t from a sample (X i,t ) 1≤i≤N , associated with the latent variables (Z i,t ) 1≤i≤N and the normalised IS weights (ω i,t ) 1≤i≤N defined in (6), we denote by (α t+1,N , θ t+1,N ) the updated value of the mixture parameters.

To approximate ( 10) and ( 11), Douc et al. (2007a) proposed the following update rule:

α t+1,N d = N i=1 ωi,t 1{Z i,t = d} , θ t+1,N d = arg max θ d N i=1 ωi,t 1{Z i,t = d} log q d X i,t ; θ t,N d . ( 12
)
The computational cost of this update is of order N whatever the number D of components is, since the weight and the parameter of each component are updated based only on the points that were actually generated from this component. However, this observation also suggests that (12) may be highly variable when N is small and/or D becomes larger. To make the update more robust, we here propose a simple Rao-Blackwellisation step that consists in replacing 1{Z i,t = d} with its conditional expectation given

X i,t , that is, ρ d X i,t ; α t,N , θ t,N
defined in (7). The resulting parameters update is given by ( 8), which we selected for Algorithm 1. Examining (7) indicates that the evaluation of the posterior probabilities ρ d (X i,t ; α t,N , θ t,N ) does not represent a significant additional computation cost, given that the denominator of this expression has already been computed when evaluating the IS weights according to (6). The most significant difference between (8) and ( 12) is that, with the former, all points contribute to the updating of the d-th component, for an overall cost proportional to D × N . Note however that in many applications of interest, the most significant computational cost is associated with the evaluation of π-which is performed exactly N times per iteration-so that the cost of the update is mostly negligible, even with the Rao-Blackwellised version.

Convergence of the M-PMC Algorithm

Convergence of the estimated parameters as N increases can be established using the same approach as in Douc et al. (2007a,b), relying mainly on the convergence property of triangular arrays of random variables (see Theorem A.1 in Douc et al., 2007a). For the Rao-Blackwellised version, assuming that for all θ's, π(q d (•; θ d ) = 0) = 0, for all α's and θ's, ρ d (•; α, θ) log q d (•, θ d ) ∈ L 1 (π), and some (uniform in x) regularity conditions on q d (x; θ) viewed as a function of θ, yield

α t+1,N d P → α t+1 d , θ t+1,N d P → θ t+1 d
when N goes to infinity. Note that we do not expand on the regularity conditions imposed on q d since, for the algorithm to be efficient, we definitely need a closed-form expression on the parameter updates. It is then easier to deal with the convergence of the approximation of these update formulas on a case-by-case basis, as will be seen in the Gaussian example of Section 3. As a practical criterion for monitoring the convergence of the algorithm we recommend computing the normalised perplexity exp(H t,N )/N (see Algorithm 1) and to interrupt adaptation when it stabilises and/or becomes sufficiently close to 1. Note that in referring to exp(H t,N ) (exponential of the Shannon entropy expressed in nat) as the perplexity, we follow the terminology in use in the field of natural language processing. The connection between the perplexity and the entropy criterion ( 2) is revealed by writing

exp -E(π, q (α,θ) ) = exp -log π unn (x) q (α,θ) (x) π(x)dx π unn (x)dx , (13)
where π unn refers to the unnormalised version of π which is effectively computable. Estimating the first integral in ( 13) by self-normalised IS as

- N i=1 ωi,t log π unn (X i,t ) q (α t,N ,θ t,N ) (X i,t )
and the second one by classical IS, as

1/N N i=1 π unn (X i,t )/q (α t,N ,θ t,N ) (X i,t ),
indeed shows that exp(H t,N )/N is a consistent estimator of exp[-E(π, q (α t,N ,θ t,N ) )]. The entropy of the IS weights is frequently used as a criterion for assessing the quality of an IS sample-together with the so-called Effective Sample Size (ESS) [START_REF] Chen | Predictive updating method and Bayesian classification[END_REF][START_REF] Doucet | Sequential Monte Carlo Methods in Practice[END_REF][START_REF] Cappé | Inference in Hidden Markov Models[END_REF]. To the best of our knowledge, however the strong connection between this criterion and the performance measure E(π, q (α t,N ,θ t,N ) ) used in the present work had not been noted before.

Variance Estimation

For the sake of completeness, we recall here the formula by which it is possible to estimate, from the IS sample, the asymptotic variance of the IS estimate. If one considers a test function h of interest, the self-normalised IS estimation of its expectation under π is π(h) = N i=1 ωi h(X i ) and its asymptotic variance is given by

υ(h) = {h(x) -π(h)} 2 π 2 (x)/q α,θ (x)dx , under the assumption that (1 + h 2 (x))π 2 (x)/q α,θ (x)dx < ∞.
The asymptotic variance υ(h) may thus be consistently estimated by [START_REF] Geweke | Bayesian inference in econometric models using Monte Carlo integration[END_REF]. Obviously, for a given function h, there is no direct link between υ(h) and the entropy criterion in (4). However it is easily shown that sup

N N i=1 ω2 i {h(X i ) -π(h)} 2 (
{h: h-π(h) ∞≤M } υ(h) = M 2 π 2 (x)/q (α,θ) (x)dx,
where the latter integral term is lower and upper bounded by 1 and exp E(π, q (α,θ) ) respectively, by direct applications of Jensen's inequality. Hence minimising E(π, q (α,θ) ) indeed reduces the worst case performance of the IS approach, at least for bounded functions (see [START_REF] Douc | Minimum variance importance sampling via population Monte Carlo[END_REF] for further discussions of this issue).

3 The Gaussian mixture case As a first example, we consider the case of p-dimensional Gaussian mixture IS densities of the form

q d (X; θ d ) = {(2π) p |Σ d |} -1/2 exp - 1 2 (X -µ d ) T Σ -1 d (X -µ d ) ,
where θ d = (µ d , Σ d ) denotes the parameters of the d-th Gaussian component density. This parametrisation of the IS density provides a general framework for approximating multivariate targets π and the corresponding algorithm is a straightforward instance of the general framework discussed in the previous section.

Update formulas

The integrated update formulas are obtained as the solution of

θ t+1,N d = arg min θ E X π ρ d (X; α t , θ t ) log |Σ d | + (X -µ d ) T Σ -1 d (X -µ d ) .
It is straightforward to check that the infimum is reached when, for d ∈ {1, . . . , D},

µ t+1 d = E X π [ρ d (X; α t , θ t )X] E X π [ρ d (X; α t , θ t )]
,

and

Σ t+1 d = E X π ρ d (X; α t , θ t )(X -µ t+1 d )(X -µ t+1 d ) T E X π [ρ d (X; α t , θ t )]
.

At iteration t of the M-PMC algorithm, both the numerator and the denominator of each of the above expressions are approximated using self-normalised importance sampling. Denoting 1{Z i,t = d} by ξ i,t , the following empirical update equations are obtained for the basic updating strategy (12):

α t+1,N d = N i=1 ωi,t ξ i,t , µ t+1,N d = N i=1 ωi,t ξ i,t X i,t X i,t N i=1 ωi,t ξ i,t = N i=1 ωi,t ξ i,t X i,t α t+1,N d , Σ t+1,N d = N i=1 ωi,t ξ i,t (X i,t -µ t+1,N d )(X i,t -µ t+1,N d ) T α t+1,N d . ( 14
)
For the Rao-Blackwellised scheme of Algorithm 1, the update is formally identical to the one above upon replacing ξ i,t by its conditional expectation

ξ RB i,t = ρ d (X i,t ; α t,N , θ t,N ) . (15)
Note that, as discussed in Section 2.2, establishing the convergence of the parameter update in this Gaussian setting will only require the assumption that ρ d (x; α, θ)x 2 is integrable with respect to π (see Douc et al., 2007a).

A simulation experiment

To illustrate the results of the algorithm presented above, we consider a toy example in which the target density consists of a mixture of two multivariate Gaussian densities. The appeal of this example is that it is sufficiently simple to allow for an explicit characterisation of the attractive points for the adaptive procedure, while still illustrating the variety of situations found in more realistic applications. In particular, the model contains an attractive point that does not correspond to the global minimum of the entropy criterion as well as some regions of attraction that can eventually lead to a failure of the algorithm. The results obtained on this example also illustrate the improvement brought by the Rao-Blackwellised update formulas in (15). The target π is a mixture of two p-dimensional Gaussian densities such that π(x) = 0.5N (x; -su p , I p ) + 0.5N (x; su p , I p ) , when u p is the p-dimensional vector whose coordinates are equal to 1 and I p stands for the identity matrix.

In the sequel, we focus on the case where p = 10 and s = 2. Note that one should not be misled by the image given by the marginal densities of π: in the ten dimensional space, the two components of π are indeed very far from one another. It is for instance straightforward to check that the Kullback-Leibler divergence between the two components of π, D { N (su p , I p ) N (-su p , I p )}, is equal to 1 2 2su p 2 = 2s 2 p, that is 80 in the case under consideration. In particular, were we to use one of the components of the mixture as an IS density for the other, we know from the arguments exposed at the end of Section 2 that the normalised perplexity of the weights would eventually tend to exp(-80). This number is so small that, for any feasible sample size, using one of the component densities of π as an IS instrumental density for the other component or even for π itself can only provide useless biased estimates.

The initial IS density q 0 is chosen here as the isotropic ten-dimensional Gaussian density with a covariance matrix of 5I p . The performances of q 0 as an importance sampling density, when compared to various other alternatives, are fully detailed in Table 1 below but the general comment is that it corresponds to a poor initial guess which would provide highly variable results when used with any sample size under 50, 000.

Proposal

N-PERP N-ESS σ 2 (x 1 )

q 0 † 6.5E
(x) = x 1 .
Quantities marked with a dagger sign are straightforward to determine, all others have been obtained using IS with a sample size of one million.

In addition to figures related to the initial IS density q 0 , Table 1 also reports performance obtained with the best fitting Gaussian IS density (with respect to the entropy criterion), which is straightforwardly obtained as the centred Gaussian density whose covariance matrix matches the one of π, that is, I p + s 2 u p u T p . Of course the best possible performance achievable with a mixture of two Gaussian densities, always with the entropy criterion, is obtained when using π as an IS density (second line of Table 1). Finally both final lines of Table 1 report the best fit obtained with IS densities of the form 0.9 D d=1 α d N (µ d , Σ d ) + 0.1q 0 (•) when, respectively, D = 1 and D = 2 (further comments on the use of these are given below). As a general comment on Table 1, note that the variations of the perplexity of the IS weights, of the ESS and of the asymptotic variance of the IS estimate for the coordinate projection function are very correlated. This is a phenomenon that we have observed on many examples and which justifies our postulate that minimising the entropy criterion does provide very significant variance reductions for the IS estimate of "typical" functions of interest.

In this example, one may categorise the possible outcomes of adaptive IS algorithms based on mixtures of Gaussian IS densities into mostly four situations: Disastrous (D.) After T iterations of the M-PMC scheme, q (α T ,θ T ) is not a valid IS density (in the sense that the importance sampling unbiasedness property does not hold due to support restrictions) and may lead to inconsistent estimates. Typically, this may happen if q (α T ,θ T ) becomes much too peaky with light tails.

As discussed above, it will also practically be the case if the algorithm only succeeds in fitting q (α T ,θ T ) to one of both Gaussian modes of π. Another disastrous outcome is when the direct application of the adaptation rules described above leads to numerical problems, usually due to the poor conditioning of some of the covariance matrices Σ d . Rather than fixing these issues by ad-hoc solutions (eg. diagonal loading), which could nonetheless be useful in practical applications, we consider below more principled ways of making the algorithm more resistant to such failures.

Mediocre (M.) After adaptation, q (α T ,θ T ) is not significantly better than q 0 in terms of the performance criteria displayed in Table 1 and, in this case, the adaptation is useless.

Good (G.) After T iterations, q (α T ,θ T ) selects the best fitting Gaussian approximation (second line of Table 1) which already provides a very substantial improvement as it results in variance reductions by about four orders of magnitude for typical functions of interest.

Excellent (E.) After T iterations, q (α T ,θ T ) selects the best fitting mixture of two Gaussian densities, which in this somewhat artificial example corresponds to a perfect fit of π. Note, however that, the actual gain over the previous outcome is rather moderate with a reduction of variance by a factor less than four.

Of course, a very important parameter here is the IS sample size N : for a given initial IS density q 0 , if N is too small, any method based on IS is bound to fail, conversely when N gets large all reasonable algorithms are expected to reach either the G. or E. result. Note that with local adaptive rules such as the ones proposed in this paper, it is not possible to guarantee that only the E. outcome will be achieved as the best fitting Gaussian IS density is indeed a stationary point (and in fact a local minimum) of the entropy criterion. So, depending on the initialisation, there always is a non zero probability that the algorithm converges to the G. situation only.

To focus on situations where algorithmic robustness is an issue, we purposely chose to select a rather small IS sample size of N = 5, 000 points. As discussed above, direct IS estimates using q 0 as IS density would be mostly useless with such a modest sample size. We evaluated four algorithmic versions of the M-PMC algorithm. The first, Plain M-PMC, uses the parameter update formulas in ( 14) and q 0 is only used as an initialisation value, which is common to all D components of the mixture (which also initially have equal weights). Only the means of the components are slightly perturbed to make it possible for the adaptation procedure to actually provide distinct mixture components. One drawback of the plain M-PMC approach is that we do not ensure during the course of the algorithm that the adapted mixture IS density remains appropriate for IS approximations, in particular that it provides reliable estimates of the parameter update formulas. To guarantee that the IS weights stay well behaved, we consider a version of the M-PMC algorithm in which the IS density is of the form

(1 -α 0 ) D d=1 α d N (µ d , Σ d ) + α 0 q 0
with the difference that α 0 is a fixed parameter which is not adapted. The aim of this version, which we call Defensive M-PMC in reference to the work of [START_REF] Hesterberg | Weighted average importance sampling and defensive mixture distributions[END_REF], is to guarantee that the importance function remains bounded by α -1 0 π(x)/q 0 (x), whatever happens during the adaptation, thus guaranteeing a finite variance. Since q 0 is a poor IS density, it is preferable to keep α 0 as low as possible and we used α 0 = 0.1 in all the following simulations. As detailed in both last lines of Table 1, this modification will typically slightly limit the performances achievable by the adaptation procedure, although this drawback could probably be avoided by allowing for a decrease of α 0 during the iterations of the M-PMC. The parameter update formulas for this modified mixture model are very easily deduced from ( 14) and are omitted here for the sake of conciseness. The third version we considered is termed Rao-Blackwellised M-PMC and consists in replacing the update equations ( 14) by their Rao-Blackwellised version (15). Finally, we consider a fourth option in which both the defensive mixture density and the Rao-Blackwellised update formulas are used.

All simulations were carried out using a sample size of N = 5, 000, 20 iterations of the M-PMC algorithm and Gaussian mixtures with D = 3 components. Note that we purposely avoided to chose D = 2 to avoid the very artificial "perfect fit" phenomenon. This also means that for most runs of the algorithm, at least one component will disappear (by convergence of its weight to zero) or will be duplicated, with several components sharing very similar parameters. Table 2 display the performance of the four algorithms in repeated independent adaptation runs. The most significant observation about Table 2 is the large gap in robustness between the non Rao-Blackwellised versions of the algorithm, which returned disastrous or mediocre results in about 60% of the cases, a fraction that falls bellow 20% when the Rao-Blackwellised update formulas are used. Obviously the fact that the Rao-Blackwellised updates are based on all simulated values and not just on those actually simulated from a particular mixture component is a major source of robustness of the method when the sample size N is small, given the misfit of the initial IS density q 0 . The same remark also applies when the M-PMC algorithm is to be implemented with a large number D of components. The role of the defensive mixture component is more modest although it does improve the performance of both versions of the algorithm (non Rao-Blackwellised and Rao-Blackwellised altogether), at the price of a slight reduction of the frequency of the "Excellent" outcome. Also notice that the results obtained when the defensive mixture component is used are slightly beyond those of the unconstrained adaptation (see Table 1). The frequency of the perfect or "Excellent" match is about 10% for all methods but this is a consequence of the local nature of the adaptation rule as well as of the choice of the initialisation of the algorithm. It should be stressed however that as we are not interested in modelling π by a mixture but rather that we are seeking good IS densities, the solutions obtained in the G. or E. situations are only mildly different in this respect (see Table 1). As a final comment, recall that the results presented above have been obtained with a fairly small sample size of N = 5, 000. Increasing N quickly reduces the failure rate of all algorithms: for N = 20, 000 for instance, the failure rate of the plain M-PMC algorithm drops to 7/100 while the Rao-Blackwellised versions achieve either the G. or E. result (and mostly the G. one, given the chosen initialisation) for all runs.

Robustification via mixtures of multivariate t's

We now consider the setting of a proposal composed of a mixture of p-dimensional t distributions,

D d=1 α d T (ν d , µ d , Σ d ) . ( 16
)
We here follow the recommendations of [START_REF] West | Modelling with mixtures[END_REF] and [START_REF] Oh | Integration of multimodal functions by Monte Carlo importance sampling[END_REF] who proposed using mixtures of t distributions in importance sampling. The t mixture is preferable to a normal mixture because of its heavier tails that can capture a wider range of non-Gaussian targets with a smaller number of components. This alternative setting is more challenging however and one must take advantage of the missing variable representation of the t distribution itself to achieve a closed-form updating of the parameters (µ d , Σ d ) d approximating ( 11), since a true closed-form cannot be derived.

The latent-data framework

Using the classical normal/chi-squared decomposition of the t distribution, a joint distribution associated with the t mixture proposal ( 16) is

f (x, y, z) ∝ α z |Σ z | -1/2 exp -(x -µ z ) T Σ -1 z (x -µ z )y/2ν z y (νz+p)/2-1 e -y/2 ∝ α z ϕ(x; µ z , ν z Σ z /y) ς(y; ν z /2, 1/2) ,
where, as above, x corresponds to the observable in ( 16), z corresponds to the mixture indicator, and y corresponds to the χ 2 ν completion. The normal density is denoted by ϕ and the gamma density by ς. Both y and z correspond to latent variables in that the integral of the above in (y, z) returns ( 16).

In the associated M-PMC algorithm, we only update the expectations and the covariance structures of the t distributions and not the number of degrees of freedom, given that there is no closed-form solution for the later. In that case, θ d = (µ d , Σ d ) and, for each d = 1, . . . , D, the number of degrees of freedom ν d is fixed. At iteration t, the integrated EM update of the parameter will involve the following "E" function

Q{(α t , θ t ), (α, θ)} = E X π E Y,Z (α t ,θ t ) { log(α Z ) + log(ϕ(X; µ Z , ν Z Σ z /Y ))| X} ,
since the χ 2 part does not involve the parameter θ = (µ, Σ). Given that

Y, Z|X, θ ∼ f (y, z|x) ∝ α z ϕ(x; µ z , ν z Σ z /y) ς(y; ν z /2, 1/2) ,
we have that

Y |X, Z = d, θ ∼ Ga (ν d + p)/2, 1 2 1 + (X -µ d ) T Σ -1 d (X -µ d )/ν d and therefore Q{(α t , θ t ), (α, θ)} = E X π D d=1 ρ d (X; α t , θ t ) log(α d ) - 1 2 E X π D d=1 ρ d (X; α t , θ t ) log |Σ d | + (X -µ d ) T Σ d -1 (X -µ d ) × ν d + p ν d + (X -µ t d ) T (Σ t d ) -1 (X -µ t d )
,

where we have used both the notation,

ρ d (X; α t , θ t ) = P α t ,θ t (Z = d|X) = α t d t(x; ν d , µ t d , Σ t d ) D =1 α t t(x; ν , µ t , Σ t )
, with t(x; ν, µ, Σ) denoting the T (ν, µ, Σ) density, and the fact that

γ d (X; θ t ) = E Y θ t {Y /ν d |X, Z = d} = ν d + p ν d + (X -µ t d ) T (Σ t d ) -1 (X -µ t d )
.

Therefore, the "M" step of the integrated EM update is

α t+1 d = E X π ρ d (X; α t , θ t ) , µ t+1 d = E X π [ρ d (X; α t , θ t )γ d (X; θ t )X] E X π [ρ d (X; α t , θ t )γ d (X; θ t )] , Σ t+1 d = E X π ρ d (X; α t , θ t )γ d (X; θ t )(X -µ t+1 d )(X -µ t+1 d ) T E X π [ρ d (X; α t , θ t )]
.

While the first update is the generic weight modification (10), the latter formulae are (up to the integration with respect to X) essentially those found in [START_REF] Peel | Robust mixture modelling using the t distribution[END_REF] for a mixture of t distributions.

Parameter update

As in Section 3.1, the empirical update equations are obtained by using self-normalised IS with weights ωi,t given by ( 6) for both the numerator and the denominator of each of the above expressions. The Rao-Blackwellised approximation based on (8) yields

α t+1,N d = N i=1 ωi,t ρ d (X i,t ; α t,N , θ t,N ) , µ t+1,N d = N i=1 ωi,t ρ d (X i,t ; α t,N , θ t,N ) γ d (X i,t ; θ t,N ) X i,t N i=1 ωi,t ρ d (X i,t ; α t,N , θ t,N ) γ d (X i,t ; θ t,N ) , Σ t+1,N d = N i=1 ωi,t ρ d (X i,t ; α t,N , θ t,N ) γ d (X i,t ; θ t,N ) (X i,t -µ t+1,N d )(X i,t -µ t+1,N d ) T N i=1 ωi,t ρ d (X i,t ; α t,N , θ t,N )
, while the standard update equations, based on (12), are obtained by replacing ρ d (X i,t ; α t,N , θ t,N ) by 1{X i,t = d} in the above equations.

Pima Indian example

As a realistic if artificial illustration of the performances of the t mixture ( 16), we study the posterior distribution of the parameters of a probit model. The corresponding dataset is borrowed from the MASS library of R (R Development Core Team, 2006). It consists in the records of 532 Pima Indian women who were tested by the U.S. National Institute of Diabetes and Digestive and Kidney Diseases for diabetes. Four quantitative covariates were recorded, along with the presence or absence of diabetes. The corresponding probit model analyses the presence of diabetes, i.e.

P β (y = 1|x) = 1 -P β (y = 0|x) = Φ(β 0 + x T (β 1 , β 2 , β 3 , β 4 ))
with β = (β 0 , . . . , β 4 ), x made of four covariates, the number of pregnancies, the plasma glucose concentration, the body mass index weight in kg/(height in m) 2 , and the age, and Φ corresponds to the cumulative distribution function of the standard normal. We use the flat prior distribution π(β|X) ∝ 1; in that case, the 5-dimensional target posterior distribution is such that

π(β|y, X) ∝ 532 i=1 Φ{β 0 + (x i ) T (β 1 , β 2 , β 3 , β 4 )} yi 1 -Φ{β 0 + (x i ) T (β 1 , β 2 , β 3 , β 4 )} 1-yi
where x i is the value of the covariates for the i-th individuals and y i is the response of the i-th individuals. We first present some results for N = 10, 000 sample points and T = 500 iterations on Figures 123, based on a mixture with 4 components and with the degrees of freedom chosen as ν = (3, 6, 9, 18), respectively, when using the non Rao-Blackwellised version (12). The unrealistic value of T is chosen purposely to illustrate the lack of stability of the update strategy when not using the Rao-Blackwellised version. Indeed, as can be seen from Figure 1, which describes the evolution of the µ d 's, some components vary quite widely over iterations, but they also correspond to a rather stable overall estimate of β, N i=1 ωi,T β i,T , equal to (-5.54, 0.051, 0.019, 0.055, 0.022) over most iterations. When looking at Figure 3, the quasi-constant entropy estimate after iteration 100 or so shows that, even in this situation, there is little need to perpetuate the iterations till the 500-th.

Using a Rao-Blackwellised version of the updates shows a strong stabilisation for the updates of the parameters α d and (µ d , Σ d ), both in the number of iterations and in the range of the parameters. The approximation to the Bayes estimate is obviously very close to the above estimation (-5.63, 0.052, 0.019, 0.056, 0.022). Figures 4 and5 show the immediate stabilisation provided by the Rao-Blackwellisation step. In this example, which is quite typical in this respect, we recommend to use less than T = 10 iterations in order to reserve most of the computational effort for increasing N , which is essential during the first adaptation steps (because the initial IS density is poor) and for the accuracy of the IS approximation in the final steps of the algorithm. Comparing the plain and Rao-Blackwellised update formulas, will really depend on how costly the parameter update is-and thus on the dimension of the model-compared to the other computational costs, and in particular the evaluation of the likelihood, which mostly depends on the number of observations. In the present case, the increase in run-time due to the use the Rao-Blackwellised formulas instead of the plain ones is negligible.

Conclusions

The M-PMC algorithm provides a flexible and robust framework for adapting general importance sampling densities represented as mixtures. The extension to mixtures of t distribution broadens the scope of the method by allowing approximation of heavier tail targets. Moreover, we can extend here the remarks made in Douc et al. (2007a,b), namely that the update mechanism provides an early stabilisation of the parameters of the mixture. It is therefore unnecessary to rely on a large value of T : with large enough sample sizes N at each iteration-especially on the initial iteration that requires many points to counter-weight a potentially poor initial proposal-, it is quite uncommon to fail to spot a stabilisation of both the estimates and of the entropy criterion within a few iterations.

While this paper relies on the generic entropy criterion to update the mixture density, we want to stress that it is also possible to use a more focussed deviance criterion, namely the h-entropy E h (π, q (α,θ) ) = D(π h q (α,θ) ) , with π h (x) ∝ |h(x) -π(h)|π(x) , that is tuned to the estimation of a particular function h, as it is well-known that the optimal choice of the importance density for the self-normalised importance sampling estimator is exactly π h . Since the normalising constant in π h does not need to be known, one can derive an adaptive algorithm which resembles the method presented in this paper. It is expected that this modification will be helpful in reaching IS densities that provide a low approximation error for a specific function h, which is also a desirable feature of importance sampling in several applications.

Figure 1 :

 1 Figure 1: Pima Indians: Evolution of the components of the five µ d 's over 500 iterations plotted by pairs:(clockwise from upper left side) (1, 2), (3, 4), (4, 1) and (2, 3). The colour code is blue for µ 1 , yellow for µ 2 , brown for µ 3 and red for µ 4 . The additional dark path corresponds to the estimate of β. All µ d 's were started in the vicinity of the MLE β.

Figure 2 :

 2 Figure2: Pima Indians: Evolution of the five Σ d 's over 500 iterations plotted by pairs for the diagonal elements: (clockwise from upper left side) (1, 2), (3, 4), (4, 1) and (2, 3). The colour code is blue for Σ 1 , yellow for Σ 2 , brown for Σ 3 and red for Σ 4 . All Σ d 's were started at the covariance matrix of β produced by R glm() procedure.

Figure 3 :Figure 4 :

 34 Figure 3: Pima Indians: Evolution of the cumulated weights (top) and of the estimated entropy divergence E π [log(q α,θ (β))] (bottom).

Figure 5 :

 5 Figure 5: Pima Indians: Evolution of the cumulated weights (top) and of the estimated entropy divergence E π [log(q α,θ (β))] (bottom) for the Rao-Blackwellised version.

Table 1 :

 1

	-4	1.5E-4	37E3

Performance of various importance sampling densities in terms of N-PERP: Normalised perplexity; N-ESS: Normalised Effective Sample Size; σ 2 (x 1 ): Asymptotic variance of self-normalised IS estimator for the coordinate projection function h

Table 2 :

 2 Number of outcomes of each category for the four algorithmic versions, as recorded from 100 independent runs.

		Disastrous Mediocre Good Excellent
	Plain	55	0	33	12
	Defensive	13	51	30	6
	R.-B.	18	1	70	11
	Defensive + R.-B.	5	11	76	8

RR n°6332

Cappé & Douc & Guillin & Marin & Robert

http://www.inria.fr ISSN 0249-6399


C Vieira

F C Aguiar

Pedro Adelar Portela

J Monteiro

P Raven

N Holmes

J Cambra

N Flor Arnau

Christian Chauvin

S Loriot

; Aguiar

; Portela

; Monteiro

; Cambra

; Flor-Arnau N

; Chauvin

Minciardi Mr

; Munné

& Ferreira

Fluvial bryophytes of Mediterranean Europe: exploring the connection between core species, climatic regions and ecological determinants

HAL is

Fluvial networks are structured across a spatial hierarchy of environmental factors, where variables acting at higher spatial scales, such as climate and geological setting, are major species occurrence ecological determinants interacting with lower spatial scale factors such as discharge and water quality or geomorphologic structure. Rivers in Mediterranean-climate regions are characterized by extremely variable annual and interannual discharge regimes, presenting seasonal and predictable floods and droughts to which faunal and floral species have adapted through their evolutionary processes. Therefore, it is assumed that fluvial communities composition between Mediterranean lotic systems may exhibit great similarity. Nevertheless, this idea has never been tested for bryophyte's communities in the Mediterranean Europe since the assessment of compositional differences between communities from different fluvial regions should be achieved using biologic and environmental field data from standardized surveys, and, until the beginning of the XXI century, this dataset was not available for bryophytic composition. Therefore, we aimed to explore the compositional segregation of bryophyte communities across three Mediterranean regions using biologic and environmental field data from standardized pan-European surveys available after the implementation and intercalibration processes of the EU Water Framework Directive. TWINSPAN classification revealed four distinct groups after two divisions (fig. 3), which do not strictly separate sites according to climatic regions, indicating other underlying ecological compositional determinants. Results from both ANOSIM and SIMPER analysis (table 1 and2) validated the groups created by TWINSPAN, as well as the indicator taxa of each group.

COMMUNITY DISTRIBUTION MODELLING

III. To better understand the distribution of the TWINSPAN groups, we used a "assemble first, predict later" community level modelling approach. We used the site groups obtained in the TWINSPAN analysis and Species Distribution Modelling techniques to model the distribution of the main Mediterranean fluvial bryophyte assemblages found.

We modelled the distribution of the groups using biomod2 ensemble forecasting package in R environment. The environmental predictors used in the modelling procedure were selected from the 10 most important variables influencing bryophyte assemblages identified in the CCA analysis.

Top environmental predictors for Group A: pH; BIO5 -max temperature of warmest month; BIO9 -mean temperature of driest quarter; BIO 15 -precipitation seasonality; Mnmanganese concentration; BIO14 -Precipitation of Driest Month.

 II. To study the influence of enviromental factors on bryophyte assemblages we performed a canonical correspondence analysis (CCA; CANOCO 5) with forward selection of environmental variables to explore the influence of environmental factors on bryophyte assemblages.

 Each TWINSPAN group corresponds to a distinct community dominated by a different assemblage of bryophyte species, presenting different and recognizable ecological preferences.

 The CCA ordination supports the segregation of TWINSPAN groups and the heterogeneity of the climatic groups (Figure 4). The first axis reflects a climatic and stream water pH gradient, represented by bryophyte communities from more alkaline streams in drier conditions to those characterizing more acidic watercourses with greater rainfall. The second axis is mainly related with an altitudinal gradient (fig. 4). We utilized the following modelling algorithms: Random Forest, Artificial Neural Networks Multiple Adaptive Regression Splines, MAXENT and Classification Tree Analysis.

Only models with an AUC evaluation score ≥ 0.7 were included in the ensemble forecast.

The ensemble predictions were obtained using the weighted mean of probabilities, which means that models with greater AUC scores had more importance in the ensemble.

The maps shown correspond to modelled probability of occurrence of each group for Mediterranean Europe.

Top environmental predictors for Group B: BIO5 -max temperature of warmest month; BIO11 -mean temperature of coldest quarter; Mn -manganese concentration; Cacalcium concentration; ELEV -elevation.

Top environmental predictors for Group C: BIO5 -max temperature of warmest month; Ca -calcium concentration; Mn -manganese concentration; ELEV -elevation; BIO11 -mean temperature of coldest quarter.

Top environmental predictors for Group D: BIO11 -mean temperature of coldest quarter; BIO5 -max temperature of warmest month; ELEV -elevation; BIO9 -mean temperature of driest quarter; Ca -calcium concentration.

MAIN CONCLUSIONS

Environmental factors Community distribution modelling Bryophyte assemblages

• Even at spatial scales as big as the Mediterranean stream water chemistry, namely pH, plays a key role in shaping the distribution of fluvial bryophytes.

• The climate also shapes the distribution of these communities, in this case the key factor is the maximum temperature in the driest and/or warmest months.

• The spatial distribution obtained for each community reflect the combined spatial patterns of water chemistry and climate.

• Low predicted probability of occurrence in sites with a recorded presence can be related with differential responses to environmental factors of the species that constitute each community.

• Bryophytic composition reveals that sites of different mediterranean climatic zones share hydrologic and geomorphologic features, creating niche for shared core taxa and other regionspecific taxa responsible for each region singularity.

  dedicate this poster to Nigel Holmes, author and colleague and an authority on river plants and the biggest influence on river conservation and restoration in the UK. By demonstrating the importance of macrophytes as ecological indicators, he transformed the way that rivers were assessed both in Britain and across other parts of Europe. He will be sorely missed and fondly remembered by his many colleagues and friends for his inspiration, boundless energy, compassion, kindness, humour and youthful enthusiasm. He was in every sense a colossus who shaped modern-day river management. C. Vieira is funded by the Fundação para a Ciência e Tecnologia (FCT) under a Postdoctoral fellowship (SFRH/BPD/63741/2009) co-funded by the Programa Operacional Ciência e Inovação -2010 and Fundo Social Europeu.

  I. First, we used floristic data collected for intercalibration purposes under the European Union Water Framework Directive and data from River Habitat Survey of 474 Mediterranean stream sites in Mediterranean Europe Climatic Areas, namely: (i) 99 sites from Mediterranean South, (ii) 246 sites from Mediterranean North, and (iii) 129 sites from Mediterranean Mountains (fig. 1 and fig. 2). In order to compare bryophyte assemblages between Mediterranean climatic regions and identify indicator (core) taxa we performed: a Two-way Indicator Species Analysis (TWINSPAN classification; CAP 3.1), a One-way Analysis of Similarity (ANOSIM) between groups obtained in the TWINSPAN classification, and a Similarity Percentage Analyses (SIMPER; PRIMER 6).

Figure 3 .

 3 Figure 3. Summary dendrogram derived from TWINSPAN classification of the 474 sampling sites in the studied Mediterranean climatic areas. For each dichotomy, the bryophyte taxa used to divide the TWINSPAN groups are indicated.

Figure 4 .

 4 Figure 4. Canonical Correspondence Analysis (CCA) ordination plot representing the Mediterranean stream sites studied and environmental variables (ELEV -elevation; Ca -calcium concentration; BIO9 -mean temperature of driest quarter; BIO14 -precipitation of driest month; BIO5 -max temperature of warmest month; BIO11 -mean temperature of coldest quarter; BIO 15 -precipitation seasonality; Mn -manganese concentration; BIO3 -isothermality). Only the 10 most influential variables obtained through forward selection procedure are represented; N-Mediterranean North sites; S: Mediterranean South sites; M: Mediterranean Mountainous sites (See Figure 1);

Figure 5 Figure 1 .

 51 Figure 5. a) TWINSPAN groups to which bryophyte communities were assigned (A-D), symbolized over different mediterranean climatic areas (a) and stream water pH (b).



  Bryophyte communities A and D are usually found in the climatic area corresponding to the Mediterranean Mountains (fig 5a), in watercourses with moderate alkalinity (fig 5b), while communities B and C are more commom in the Mediterranean North and South (fig 5a), in more acidic streams (Fig. 5b).

Table 2 .

 2 Results of SIMPER analysis for bryophyte species composition of stream sites within each TWINSPAN group (A-D) (Asim = average similarity of sampling sites in the group; Av.abund = average abundance; Co% = percentage of the contribution of species for similarity of the group).

	Group A (Asim: 20.95)

Group C (Asim: 22.07) Group D (Asim: 21.05) Species Av.abund Co (%) Species Av.abund Co (%) Fissidens crassipes subsp. warnstorfii 0.44 18.33


* Biological macrophytes index for rivers (Haury et al., 2006): French macrophytes index for ecological status assessment of streams METHODOLOGY Relative abundance of bryophytes Relative abundance of bryophytes Morphodynamical gradient (32%) Size gradient (24%)

Developing macrophytes-based methods for assessing and managing the ecological status in rivers needs to deal with different analysis scales: floristic data are recorded at the site scale while the anthropic pressures are defined at the large scale of catchment area. Focusing on hydromorphology highlights the needs of multi-scales approach. Indeed these environmental parameters cover functioning (hydrology) and features (substrate, river shape) at nested scales. To be able to understand the relationship between global descriptors (anthropic pressures) and precise population composition and structuration (bioindicators), focusing the right scale for each approach is essential.

The search for sensitive macrophytes metrics for assessing the impact of environmental parameters were tackled by a PLS modeling distinguishing pressures at catchment, reach and site scales so as to respect the hierarchical spatial organization. This multi-scalar method was developed to identify and quantify the biotic response through these different scales.

Focusing on hydromorphological pressures and the local scale, a PCA analysis was realised on a biological metric in order to refine the PLS model reading. This study was undertaken using the French river monitoring databases (chemical, hydromorphological and biological parameters). The habitat description (local scale) was derived from the macrophyte survey protocol.

Correlations between PLS predictors are shown by PCA LOCAL SCALE

At the global scale the PLS shows that nutrients load and organic matter had a predominant effect on the biological compartment, and that land use features played an integrating role for the pressures expression. Oppositely, the analysis shows the difficulty to characterise the role of hydromorphological descriptors assessed at the medium scale. That is clearly shown using a global metric as IBMR* (Index value and Ecological Quality Ratio).

If at a large scale, the descriptors of pressures cannot directly explain the response of the biological metrics, the link is more clearly evidenced by analysing the local hydromorphological parameters. The habitat descriptors are recorded in the same time and at the same site where the floristic data is surveyed (100 m stretch).

At this consistent scale, data describing both population and habitat are closely related. The exemple of the bryophytes relative abundance shows a clear gradient of response of the biological metrics regarding the morphodynamical (MD) and size (S) gradient of the stream.

The Water Framework Directive required to fit the biological assessment methods as sensitive indicators of catchment-wide scale pressures. The whole results of the study shows that the direct relationship is basically poor at the site scale, as it can be expected from the ecological nested concepts. To be able to fit efficient and relevant assessment methods in the one hand and diagnosis methods in the other, developing methodologies must be undertaken at the consistent and right scale, both for biological and environmental parameters. If necessary, relationships between catchment-wide pressure descriptors and local habitats features (morphology, substrate, chemicals) must be defined in an other approach for establish proper links on the global population response model (e.g. DPSIR models).

  Macrophytes response to environmental stressors:changing scale from habitat characteristics to catchment-wide pressures S. Loriot, T. Feret, Bertrand Villeneuve, Christian Chauvin To cite this version: S. Loriot, T. Feret, Bertrand Villeneuve, Christian Chauvin. Macrophytes response to environmental stressors: changing scale from habitat characteristics to catchment-wide pressures. 14th International Symposium of Aquatic Plants, Sep 2015, Edinburgh, United Kingdom. pp.1, 2015. hal-02601867Macrophytes response to environmental stressors : changing scale from habitat features to catchment-wide pressures Sandrine LORIOT 1 , Thibaut FERET 1 , Bertrand VILLENEUVE 2 , Christian CHAUVIN 11 Irstea Bordeaux, CARMA, Aquatic Ecosystems and Global Changes, 50 avenue de Verdun -33612 Cestas, France.(contact: sandrine.loriot@irstea.fr) 2 Irstea Lyon, LHQ Aquatic Habitats, Ecology and Pollutions, 5 rue de la Doua -69626 Villeurbanne, France.Multi-scale approach : from catchment structure to site habitat 14 th International Symposium on Aquatic Plants Edinburgh

  is performed in the framework of pluri-annual project supported by ONEMA and the French Environment Ministry, targeting the development and the implementation of the WFD water ecological assessment methods.


est destinée au dépôt

Context : Buffer strips are a way to limit pesticides transfer from fields to aquatic environment. Contaminants transfers are strongly related to water pathways. Modelling water behaviour in areas may help to better understand pesticides lateral transfers and interactions between surface and subsurface flows. Saturated conductivity is one of the most influent hydrodynamic characteristics for water flows and is difficult to measure on the field. This study aims at better understanding the effect of saturated hydraulic conductivity heterogeneity on water flows with the help of a physically based model. It is based on data collected on an experimental buffer strip ; such devices are particularly active areas for what concerns water flow, both in surface and subsurface.tance Objectives : (1) Compare different ways to describe conductivity heterogeneity using geostatistics and exclude non consistant ones (mathematical criteria).

(2) Test the chosen descriptions against field observations on an artificial event under specific criteria. -the law (in this case lognormal and Gaussian) -the correlation length (8 m and 20 m) -the enforcement or not with ponctual measured conductivity points (16 points distributed over the area from surface to 40-cm depth).

Results

Conclusion & perspectives

Exemple of case F statistical fiels : with enforcement points on horizons 1 and 2. Definition of the conductivity distribution characteristics of the 7 simulated cases. A : Boxplots of maximum area affected by runoff (left) and B : Boxplot of cumulated area reached by surface water infiltration at 50 cm-depth and t = 3 hours (right) in the seven simulation cases.

Comparison of the seven cases on two outputs : one for the surface and one for the subsurface  Cases A, E and D seem to be the most consistent cases for the chosen criteria.

Comparison of cases A, E, D and G with field datas

Introduction & Objectives

Experimentation field (Saint-Joseph buffer strip, Beaujolais) with piezometers positions.

Experimentation field (6 m-long and 4 m-wide).

Mean values are stable for all distributions

For no enforcement but with the correlation length cases, results are very variable both for surface and subsurface.

The gaussian case results have a low variance for surface but are very fluctuating for subsurface.

With an homogeneous soil, results are close to the mean of other cases.

3

Additional information : microtopography correlated with first horizon conductivity.

Watertable dynamics at piezometer 2.

Homogeneous case is close to observation data.

With enforcement and correlation length cases (E and F), results are scattered but globally follow the observed dynamic.

Case A without correlation is closer to observed data than cases E and F. Cases A watertable behaviour is comparable with homogeneous case.

 With a heteregeneous distribution, taking into account conductivity measured points improves results stability.  At this scale, an horizontal length distribution of 8 m or 20 m does not stabilize the results by creating low or high conductivity zones.  In this precise context, the two distributions that lead to the closest results to observed data are the less heterogeneous ones : lognormal law without length correlation and homogeneous ones. Further perspective : test those distributions on a natural event to verify the conclusions and eventually change the scale of study (distributions with correlation length could benefit from a larger scale (Binley 1989)).

  Influence of soil hydrodynamic characteristics variability on surface and subsurface flows at a vegetative buffer strip scale Laura Gatel, Claire Lauvernet, Claudio Paniconi, Nadia Carluer, Etienne Leblois To cite this version: Laura Gatel, Claire Lauvernet, Claudio Paniconi, Nadia Carluer, Etienne Leblois. Influence of soil hydrodynamic characteristics variability on surface and subsurface flows at a vegetative buffer strip scale. IRTG "Integrated Hydrosystem Modelling", Apr 2015, Tübingen, Germany. pp.1, 2015. hal-Influence of soil hydrodynamic characteristics variability on surface and subsurface flows at a vegetative buffer strip scale. Laura Gatel a* , Claire Lauvernet a , Claudio Paniconi b , Nadia Carluer a , Etienne Leblois a . a-Irstea, UR MALY, 5 rue de la Doua, 69100 Villeurbanne, France (*corresponding author: laura.gatel@irstea.fr), c-INRS, centre Eau-Terre-Environnement, 490 rue de la Couronne, Québec (Québec) G1K 9A9, Canada.

MaterielModel:

  CATHY (CATchment HYdrology) (INRS-ETE Québec) is a physically-based 3-D model (Camporese 2010) for water transfer in surface and subsurface (using respectively the diffusion wave and Richards equations). References : Binley A., Egly J., Beven K. (1989). A physically based model of heterogeneous hillslopes: 1. Runoff production. Water resources Research, 25(6):1219-1226. Camporese M., Paniconi C., Putti M., and Orlandini S. (2010). Surface and subsurface flow modeling with path-based runoff routing, boundary condition-based coupling, and assimilation of multisource observation data. Water Resources Research, 46(2). Lacas J-G. (2005). Processis de dissipation des produits phytosanitaires dans les zones tampons enherbées : étude expérimentale et modélisation en vur de limiter la contamination des eaux de surface. PhD thesis, Montpellier 2. Leblois E. and Creutin J-D. (2013). Space-time simulation of intermittent rainfall with prescribed advection field : Adaptation of the turning band method. Water resources Research, 49(6):3375-3387.

  Study site : 6-m long area enclosed in a 20 m-long buffer strip between a vineyard and the Morcille river (Beaujolais, France). The study will use data from an artificial experiment performed in 2006 (200 mm in 40 minutes). Geostatistical tool (Leblois & Creutin, 2013) : Generation of statistical fields of conductivity. The user can choose :

  mesh of the area as represented with Cathy with boundary conditions.


Nicolas Forquet

Matthieu Dufresne

Simple deterministic model of the hydraulic buffer effect in septic tanks

Keywords:

Septic tanks are widely used in on-site wastewater treatment systems. In addition to anaerobic pretreatment, hydraulic buffering is one of the roles attributed to septic tanks. However there is still no tool for assessing it, especially in dynamic conditions. For gravity fed-system, it could help both researchers and system designers. This technical note reports a simple mechanistic model based on the assumption of flow transition between the septic tank and the outflow pipe. The only parameter of this model was calibrated using CFD modeling for a wide range of discharge rates. The resulting model highlights that a septic tank plays a hydraulic buffer role when faced with sudden and large discharge flow but this role tends to disappear when input hydrographs are smoother. In those cases there is an observable lag between the input hydrograph and outflow hydrograph.

Introduction

On-site wastewater systems usually consist of a septic tank followed by a treatment unit. In many countries (including the USA and France), treatment units are commonly gravity-fed from the outlet of the septic tank (e.g. a drainfield trench or vertical flow sand filter). Because flow at the outlet of a septic tank is not constant, the distribution over the surface of the treatment unit is rather uneven, especially in early filter operation [START_REF] Bridson-Pateman | The effects of dosed versus gravity-fed loading methods on the performance and reliability of contour trench disposal fields used for on-site wastewater treatment[END_REF][START_REF] Gill | Nutrient loading on subsoils from on-site wastewater effluent, comparing septic tank and secondary treatment systems[END_REF]. Flow variability stems mainly from household water usage patterns, but the septic tank also induces some flow modulation. The hydraulic buffer effect of septic tanks is often cited but to our knowledge there is still no method for quantifying it. It is often assumed that for the purposes of studying of long-term phenomena (like clogging of the treatment unit at month-long or year-long scale), the outflow can be considered constant (Winstanley & Fowler, 2013). However, as we gain progressively more knowledge on the actual hydrograph produced by a household, it may be interesting to quantify how effectively the septic tank can buffer large inflow rates. This may be of particular interest for estimating the efficiency of gravity-driven distribution on secondary treatment unit and to eventually optimize it.

1 Irstea, UR MALY, centre de Lyon-Villeurbanne, 5 rue de la Doua -CS 70077, 69626 Villeurbanne Cedex, France. nicolas.forquet@irstea.fr 2 Ecole Nationale du Génie de l'Eau et de l'Environnement de Strasbourg, Laboratoire ICube (Université de Strasbourg, CNRS, INSA de Strasbourg, ENGEES), ICube, 2 rue Boussigault 67000 Strasbourg, France.

Two possible approaches were identified: (i) using a tank model where the law governing outflow is obtained by a statistical learning or neural network method (Vazquez et al., 1999) or (ii) the overflow analogy. A statistical learning or neural network method is able to mimic complex hydraulic systems without the need to compile advanced knowledge of the constitutive elements. However, it requires a large dataset for model learning, which is not compatible with our needs. The overflow analogy is based on a simple mechanistic approach based on the assumption that critical flow occurs at the outlet of the septic tank. In this paper, we briefly present the model and its practical implementation, and then report selected results based on several hydrographs.

Model presentation

Model equations

Figure 1 presents a sketch of the usual design of a septic tank. The upstream section of the outflow pipe of a septic tank is a local maximum of the bottom profile. With no downstream influence (guaranteed by the large slope of the pipe, which is usually over 0.5%; AFNOR, 2007), this configuration is responsible for producing critical flow [START_REF] Hager | OpenFOAM -The open source CFD toolbox -User guide[END_REF]. Critical flow is the transition between subcritical flow (here, a nearly horizontal water surface with a very small velocity in the septic tank) and supercritical flow (here, a fast flow in the outflow pipe). The occurrence of critical flow guarantees a direct relationship between water level h in the tank and the outflow discharge (Q out [L 3 T -1 ]). We used this relationship in association with a water mass balance in the septic tank to build a time-dependent model. The discharge Q out corresponding to the critical water depth h c [L] can be evaluated considering a Froude number equal to unity (transition from subcritical to supercritical flow), according to:

hc c out gD S Q  (1)
where

Q out is outflow discharge [L 3 T -1 ], S c is critical cross-section [L 2 ], g is gravitational acceleration [L T -2 ] and D hc is critical hydraulic diameter [L].
Both S c and D hc are linked to the critical water depth (h c ) by relationships (Equations 2, 3 and 4) using the angle δ c as illustrated in Figure 2.

Author-produced version of the article published in Water and Environment Journal, 2015, 29(3), 360-364. The original publication is available at : http://onlinelibrary.wiley.com/doi/10.1111/wej.12114/pdf doi : 10.1111/wej.12114

h c D  c
  c c c c D S    cos sin 4 2   (2)   c c c c hc D D     sin 4 cos sin   (3)   c c D h  cos 1 2   (4)
where D is outlet pipe diameter [L]. The angle δ c may be expressed as a function of D and h c :

        D h a c c 2 1 cos  (5)
Equation 1 links the outflow discharge to the critical water depth in the outflow pipe. We are now seeking out a relationship between critical water depth (h c ) and the water depth in the tank (h [L],

measured from the invert of the outflow pipe). Knowing the critical water depth (h c ), the critical energy head H c [L] can be calculated with the following expression.

2 2 2 c out c c gS Q h H   (6)
The head loss ΔH [L] between the tank and the critical section can be evaluated as a local head loss, thus:

2 2 2 c out gS Q K H   (7)
The loss coefficient K [-] was evaluated by CFD using the OpenFOAM software package (2013) for a 100 mm diameter pipe and a flow rate ranging from 0.10 to 1.50 L/s. The conclusion of the numerical simulations is that the loss coefficient K is approximately 0.4 for the whole discharge range. Using this value and an estimation of the numerical uncertainty based on a grid sensitivity analysis, the uncertainty on the outflow discharge (Q out ) for a given water depth (h) was evaluated as 5%. Finally, the water depth in the tank (h) can be evaluated using equation 7 (Bernoulli equation written between the critical section and the tank where the velocity head is close to zero).

H H h

c    (8)
Based on equations 6 and 7, equation 8 can be rewritten as:

  2 2 2 1 c out c gS Q K h h    (9)
Replacing h c and S c in equation ( 9) by expressions dependent solely on Q out will lead to an expression relating tank water depth (h) to outflow discharge (Q out ). However, analytically solving this equation would prove cumbersome due to the sinusoidal functions involved. An alternative solution was found that consisted in rewriting equation 8 into a minimization problem. Incorporating equation ( 1) into equation ( 9) implies:

hc c D K h h 2 1    (10)
Rewriting equation 10 into the form of an objective function (obj.fun) gives:

2 2 1 ) ( .           hc c c D K h h h fun obj (11)
For a given value of h, finding the value of h c that minimizes the objective function makes it possible to compute the outflow discharge (Q out ) using equation 1. Once this relation has been established, it can be associated to the septic tank mass balance equation to build a time-dependent model. For a septic tank, the water mass balance can be written as:

out in Q Q dt dh S   (12)
where S is horizontal surface of the septic tank at the level of the invert of the outflow pipe [L 2 ], and

Q in is inflow rate [L 3 T -1 ].
After an explicit discretization, equation 12 becomes:

) ( ) ( ) ( ) ( t Q t Q t t h t t h S out in       (13)
Author-produced version of the article published in Water and Environment Journal, 2015, 29(3), 360-364. The original publication is available at : http://onlinelibrary.wiley.com/doi/10.1111/wej.12114/pdf doi : 10.1111/wej.12114

The value of Q in (t) is an input while the value of Q out (t) needs to be estimated (except for the initial value). This is carried out by minimizing the objective function (11) at time t.

Model implementation

The model was implemented in R (R Core Team, 2013). The method used for minimizing the objective function is a combination of golden section search and successive parabolic interpolation [START_REF] Brent | Algorithms for Minimization without Derivatives[END_REF]. Accuracy was set to 1 x 10 -10 m.

Practical test cases

This technical note presents three test cases. [START_REF]Mise en oeuvre des dispositifs d'assainissement non collectif[END_REF] It is often assumed that the largest inflow rate for a septic tank (if properly disconnected from rainwater) corresponds to the emptying of a bath (200 litres in 3 minutes). ( 2) In France, new treatment systems require authorization before being commercialized. Since 2009, this authorization is given based on the results of normalized experiments (AFNOR, 2013) carried out by accredited laboratories. Feeding of the wastewater treatment system during these normalized experiments is carried out according to a distribution of the daily hydraulic load that is based on the assumed consumption of a typical household. Based on this distribution, a synthetic hydrograph was generated, corresponding to a two person-household. Daily hydraulic load was estimated at 84 L/pers./day according to [START_REF] Cauchi | Characteristics of raw water from the individual house[END_REF]. (3) [START_REF] Butler | Modeling dry weather wastewater flow in sewer networks[END_REF] and [START_REF] Butler | Synthesising dry weather flow input hydrographs: A Maltese case study[END_REF] presented synthetic hydrographs of wastewater discharge in person-equivalents. Despite the fact that these studies were done on sewage, they are often cited as a benchmark for on-site wastewater treatment systems due to the lack of input hydrographs in this

Author-produced version of the article published in Water and Environment Journal, 2015, 29(3), 360-364. The original publication is available at : http://onlinelibrary.wiley.com/doi/10.1111/wej.12114/pdf doi : 10.1111/wej.12114 domain [START_REF] Roland | Comparative analyses of seepage: clogging tools for the diagnosis[END_REF]. Here, we used the one presented in [START_REF] Butler | Synthesising dry weather flow input hydrographs: A Maltese case study[END_REF] as an input in our model. The daily hydrograph was pre-normalized so that daily hydraulic load is the same as in case 2. For all three test cases, the characteristics of the septic tank are the same, i.e. a 4 m 2 area at the invert of the outflow pipe (100 mm in diameter) that corresponds to a commercial standard.

Figure 4 presents the results. A single event, like a bath emptying (case a), is considerably smoothed by the septic tank. The outflow discharge reaches only 57% of the inflow rate at the end of loading.

Many treatment systems (e.g. sand filters) downstream of the septic tank have gravity-driven distribution systems and therefore depend on inflow velocity to ensure an even distribution. If the septic tank significantly smoothes its outflow discharge, this could affect the quality of hydraulic distribution in the treatment system. When applied to more averaged hydrographs such as those obtained synthetically based on daily hydraulic load and typical household water usage (AFNOR, 2013; case b) or those experimentally observed at sewer level [START_REF] Butler | Synthesising dry weather flow input hydrographs: A Maltese case study[END_REF]; case c), the smoothing effect tends to be less significant: mean outflow discharges are only 0.6% and 3.75% lower than mean inflow rates for the second and third test cases, respectively. In these cases, the septic tank only induces a lag in the propagation of the hydrograph. In addition, case b shows that the septic tank may help ensure a better spreading of hydraulic load over the day.

Evaluation of septic tank buffering on effluent distribution over the secondary treatment unit would requires better measurements of the inlet hydrograph with a small time resolution. The hydrograph suggested by (AFNOR, 2013) is too averaged and actually presents a shape close to the observed in sewers. [START_REF] Patel | A comparison of gravity distribution devices used in onsite domestic wastewater treatment systems[END_REF] measured outflow at the inlet of gravity distribution devices (after the septic tank) and concluded that the most common flow rates were between 0.0016-2.0.03 l h -1 with peak values up to 0.2 l h -1 . For secondary treatments that are not fed by gravity, the septic tank buffer effect may help to ensure a better spreading of the influent over time. Furthermore, the current model could be easily applied to septic tank alternative designs such as tank in series and tank in parallels as they are already widespread in USA and tend to develop in Europe. A sensitivity analysis on parameters S and D could also indicate which one influence the most buffer effect and lag time.

Finally, we would like to stress that not only septic tank induces a buffer effect. Pipes, conducting flow to and out of the septic tank, may also be of importance regarding their diameter (typically 100 mm).

Conclusions

We built a simple mechanistic model suitable for modeling the hydraulic buffering induced by a septic tank. The only constant in the model, i.e. the loss coefficient, was calibrated using deterministic CFD modeling. Results highlight significant modulation of flow in single relatively large events such as a bath emptying. However, on smoother hydrographs, such as those currently available for characterizing a household effluent, the flow modulation only adds a lag time into the influent hydrograph. As we progressively gain more knowledge on the actual shape and amplitude of the hydrograph at the inlet of a septic tank, this simple tool could prove be useful for modeling septic tank outflow and its impact on the spread of wastewater over the treatment unit in configurations based on gravity-driven distribution.

Author-produced version of the article published in Water and Environment Journal, 2015, 29(3), 360-364. The original publication is available at : http://onlinelibrary.wiley.com/doi/10.1111/wej.12114/pdf doi : 10.1111/wej.12114

Figure 1 .

 1 Figure 1. Sketch of the water flow through a septic tank

Figure 2 .

 2 Figure 2. Relationships between critical water depth and hydraulic section and pipe diameter

Figure

Figure

  Figure 3. Computation steps

Figure 4 .

 4 Figure 4. Inflow and outflow hydrographs for three test cases: a) bath emptying, b) synthetic daily dataset, c) dry weather sewage flow[START_REF] Butler | Synthesising dry weather flow input hydrographs: A Maltese case study[END_REF]

Author-produced version of the article published inWater and Environment Journal, 2015, 29(3), 360-364. The original publication is available at : http://onlinelibrary.wiley.com/doi/10.1111/wej.12114/pdf doi : 10.1111/wej.12114


A Gao

Marina Dabrin

Coquery

DGT for the monitoring of trace metals, mercury and arsenic in surface freshwaters

In situ intercomparison exercise on "home-made"

00er

In situ intercomparison exercise on "Home-made" DGT for the monitoring of trace metals, mercury and arsenic in surface water A. Magnier 1 , L. Lesven 2 , Y. Gao 3 , A. Dabrin 1 , M. Coquery 1 1 Irstea, UR MALY, 5 rue de la Doua-CS70077, F-69626 Villeurbanne cedex, France 2 Laboratoire LASIR (UMR 8516), Equipe Physico-chimie de l'Environnement, Université Lille 1, F-59655 Villeneuve d'Ascq, France 3 Department of Analytical and Environmental Chemistry (ANCH), Vrije Universiteit Brussel (VUB), Brussels, Belgium Most monitoring programs for surface waters within a legislative framework such as the Water Framework Directive (WFD) rely on conventional techniques such as spot or automated sampling to determine total dissolved concentration of a target substance. These sampling methods are generally time-consuming, do not take into account temporal variability and could induce samples contamination, loss of analyte or speciation modifications. The Diffusive Gradient in Thin films (DGT) technique is an alternative sampling method to assess a time-weighted average (TWA) metal concentration in surface waters, as it takes into account metal variations during the period of exposure. This tool coupled with different resins adapted to different kind of substances is commercially available but several laboratories are now developing their home-made DGTs. As each laboratory prepares its home-made gel or resin in different working conditions by using different reagents (i.e. purity and/or supplier) and by using different methods to determine diffusive coefficient, the comparison of results obtained with these tools could be not relevant. To answer this question, the objective of our study was to compare commercial and "home-made" DGT performances for cationic metals (Cd, Pb, Cu, Zn), arsenic and total mercury during an in situ intercomparison exercise in Deûle River.

Mercury

Lab 2

Metals

Lab 3

Different kinds of chelating resins using in DGT method. In our study, Chelex 100 , 3-MP, spheron thiol and ferrihydrite chelating resins were used. Intercomparison exercise has shown a good fit between the different concentrations of trace metals and total mercury sampled by "home-made" DGT. Results have clearly indicated that diffusion coefficients have to be determined for each "home-made" DGT (chelating resin+diffusive gel+filter) in order to achieve the best suitable results. Working conditions from clean lab to polarographic laboratory (using mercury drop electrode) impact blank results

Metals Total mercury Arsenic

All the teams implicated in this project would like to acknowledge funding support from Aquaref and Onema.

For more informations: Ludovic.lesven@univ-Lille1.fr

Conclusion

DGT Conference 2015 San Sebastian, Spain

A:

  Tinytag used to record average temperature during DGT deployment (B) (C: average temperature : 19°C) DGT exposition in Deûle river (North of France) near Douai Town. Station situated near Metaeurop, former which has strongly impacted this site until 2003 by metal discharges. Major parameters measured in the sampling site from the beginning to the end of DGT exposition A B C Blank values (B) shown differences imputed to clean room conditions (lab 3 : polarographic environment) Same concentration range (C) for all results Different levels of Hg accumulation on 3-MP resin for each lab (B) due to « home-made » DGT used Calculation of Hg concentration using experimental diffusive coefficient (A) Analysis of arsenic, total mercury and cationic DGTs performed respectively by Lab 1, 2 and 3. Note that comparison is done between manufacturing of « home-made » chelating resin+diffusive gel and do not into account treatment + analysis with specific apparatus.  Same diffusive gel thickness for all DGTs : ~0.8 mm (A)  Same diffusion coefficient (Lancaster) (A) No difference between mass and concentration variations Metal concentration from home-made DGT close to commercial DGT (B) values < LQ (B) All arsenic concentrations determined from DGT method are differents (C) Need to characterize own diffusion coefficients for each « home-made » DGT used (A)


F Serveto

A Daval

F Desbiolles

M Coquery

C Miège

Materials and methods: an efficient multiresidue method Application to different types of wastewaters: concentration levels Conclusions

Keywords: Hypolipemiants, AC: Anticarcinogen, AD: Antidiabetic, AE: Antiepileptic, BD: Bronchodilator * : negative ionization (ESI -)

Sample filtration: Glass filber filter 0.7 µm Sample acidification to pH 5 Sample purification: OASIS HLB ® 500 mg, 6 mL Conditioning : 5 mL MeOH, 5 mL H 2 O Washing : 5 mL H 2 O Elution : 2*5 mL MeOH Analysis: UHPLC-ESI(+/-)-QqQ-MS-MS ® Shimadzu Nexera 2 , API 4000 ® ABSciex Column C18

Indicator pesticides and pharmaceuticals were selected according to:

-their frequency of quantification in untreated wastewaters, secondary and tertiary wastewater effluents and in receiving rivers; -their toxic effect in aquatic and terrestrial ecosystems; -their large range of physico-chemical properties: octanol-water partition coefficient (log Kow), acid dissociation constant (pKa) and chemical classes.

selected substances

including 4 priority substances of the WFD and 5 watch-list substances of the WFD daughter directive (2015). One of the main issue of the Water Framework Directive (WFD, 2000/60/CE) is to reach good ecological and chemical status of waters bodies. The list of WFD priority substances has been amended in 2013 (2013/39/UE) and 2015 (2015/495/UE) by identifying new substances posing a significant risk to, or via, the aquatic environment. In fact, these new substances, so called "emerging contaminants", including numerous pharmaceuticals and pesticides, are frequently detected in rivers and wastewaters in concentrations from ng/L up to several µg/L. Hence, it is now of crucial importance to monitor wastewater treatment plant (WWTP) effluents since they are known one of the main vector of micropollutants to surface water.  Our challenge was to develop a reliable analytical method to quantify emerging contaminants at trace level in various types of waters.

Application to WWTP influent, secondary and tertiary treatment effluents of an ozonation process and to influent and effluent of a Planted Discharge Area.

26 substances quantified at concentrations ranging from 1.5 to 68 000 ng/L. water sample

• This study showed that our multiresidue analytical method is suitable for the quantification of 41 emerging micropollutants from ng/L to 100 µg/L in various water matrices. • Further work is needed to better appreciate matrix effects on various aqueous matrices.

• This method will be adapted for biota matrices (Gammarus fossarum) and passive samplers extracts (Polar Organic Chemical Integrative Sampler).

  at concentrations ranging from 0.3 to 1 800 ng/L. IPU, CYCLOP and APZ were never quantified. Analytical problems for NFLUOX, ACE, CLINDA, PARA, OFLO and DIM.

  from LOQ values to 200 ng/L. -Recoveries evaluated on a surface water spiked at 60 µg/L under reproducibility conditions (n=7): • Only qualitative analysis for NORFLO, CIPRO, METFOR, AZI, FENO. • Work in progress for negative ionization compounds. Identification/quantification: MRM mode acquisition Confirmation criteria are described in the EU council decision Internal calibration Use of 22 labelled molecular analogues as internal standards and surrogates

Acknowledgments:

  We thank colleagues for sample collection; and the French National Agency for Water and Aquatic Environments (ONEMA) and the Rhône Méditerranée Corse Water Agency through the MICROPOLIS and ZRV projects for financial support. Extract dilution according to the matrix and expected concentrations A multiresidue method to evaluate emerging micropollutants levels in waters 19 substances quantified at concentrations ranging from 0.5 to 390 ng/L. DCA, DIA, DEA, SMZ, ATZ, DIAZ, APZ and DIM were never quantified. Analytical problems (matrix effects or QC out of specifications) for PIRI, PARA, IPU, ATE and ERY. of quantification; N.D.: not determined UHPLC-MS-MS system Quality controls (QC): -Laboratory blanks (source water) to check non contamination. -Spiked samples to check recoveries. -QC standard solutions at 2 ng/L and 40 ng/L every 6 samples to check accuracy of the analysis and non deviation of instrumental measure. Acceptation criteria for QC standard solution: between 70 and 130 %.


est destinée au dépôt

Complementarity of passive samplers and biomonitoring tools to evaluate surface water contamination by pesticides in agricultural catchments

COMPLEMENTARITY OF PASSIVE SAMPLERS AND BIO-INDICATORS TO EVALUATE SURFACE WATER CONTAMINATION BY PESTICIDES IN AGRICULTURAL CATCHMENTS

Le Dréau Matthieu, Chaumot Arnaud, François Adeline, Geffard Olivier, Margoum Christelle, Pesce Stephane (1) , Martin Caroline, Mazzella Nicolas (2) , Gouy Véronique (1)

(1) Irstea 5, rue de la Doua CS 70077 69626 Villeurbanne cedex -matthieu.le-dreau@irstea.fr

(2)

Irstea 50, avenue de Verdun 33612 Cestas cedex

Summary

Grab sampling may be scarcely adequate to evaluate the impact on surface waters of pesticide used in agriculture because of its poor spatiotemporal representativeness compared to the high variability of the contamination patterns. Indeed, it is still not easy to link the effect of agricultural practice changes to water quality at the catchment scale. Furthermore, biological indicators may also be useful in addition to pesticide concentration determination to evaluate their impact on ecosystems.

In this way, the contribution of new integrative tools has been evaluated. Passive samplers and biological indicators have been applied together at the catchment scale so as to identify their coherence and complementarity. The tested chemical tools are the passive-SBSE (Stir Bar Sorptive Extraction) dedicated to hydrophobic to moderately hydrophilic compounds and the POCIS (Polar Organic Chemical Integrative Sampler) dedicated to hydrophilic compounds. As concerns biological approach, it consisted in the in situ deployment of i) caging methodology with a highly ecological relevant amphipod species, Gammarus fossarum, coupled to the measurement of toxicity markers (acetylcholine esterase inhibition, feeding rate) and ii) the study of microbial and macroinvertebrate communities (measurement of leaf decomposition rates) as well as phototrophic periphytic communities (measurement of Pollution Induced Community Tolerance). These tools have been applied during june 2014 (corresponding to vineyard and corn period of application) in two contrasted agricultural catchments. The Ardières basin situated in the Beaujolais region (France) is submitted to viticulture pressure and a marked pesticide contamination has already been highlighted. Three measurement stations have been chosen so as to search for a contamination gradient from upstream to downstream. The Coise basin is located in the Monts du Lyonnais (France) in a zone of mixed farming breeding with a predominance of meadows, corn and cereals. Three measurement stations have been selected on three sub-catchments corresponding to differentiated pesticide use. In total, 49 pesticides (28 herbicides, 9 fungicides and 12 insecticides) have been analyzed as well as potential confusing parameters (metal(oïd)s, nutrients, temperature, luminosity, water flow). Results show the high coherence of chemical and biological tools responses. All of them were able to reflect the spatial contamination gradients as well as the specific contamination and impact associated to each agro-pedo-climatical context. A weekly exposure also permitted to take into account the temporal variability of the contamination according to flood occurrence. Finally, their sensitivity, reactivity and integration capacity make these tools relevant and promising to assess agricultural pesticide impact on surface water in addition to grab sampling.

  M. Le Dreau, Arnaud Chaumot, A. Francois, Olivier Geffard, C. Margoum, Stéphane Pesce, C. Martin, Nicolas Mazzella, Véronique Gouy To cite this version: M. Le Dreau, Arnaud Chaumot, A. Francois, Olivier Geffard, C. Margoum, et al.. Complementarity of passive samplers and biomonitoring tools to evaluate surface water contamination by pesticides in agricultural catchments. Symposium in Pesticide Chemistry, Sep 2015, Piacenza, Italy. pp.1, 2015. hal-02602027

Symposium in Pesticide Chemistry, Piacenza, Italy, 2-4 September, 2015

Acknowledgement

This work benefits from the technical and financial support of Onema and the Ecophyto Program as well as the ANR CESA Program (Gamma and Potomac projects).