Wojciech De Roeck

François Huveneers

Glassy dynamics in strongly anharmonic chains of oscillators

We review the mechanism for transport in strongly anharmonic chains of oscillators near the atomic limit where all oscillators are decoupled. In this regime, the motion of most oscillators remains close to integrable, i.e. quasi-periodic, on very long time scales, while a few chaotic spots move very slowly and redistribute the energy across the system. The material acquires several characteristic properties of dynamical glasses: intermittency, jamming and a drastic reduction of the mobility as a function of the thermodynamical parameters. We consider both classical and quantum systems, though with more emphasis on the former, and we discuss also the connections with quenched disordered systems, which display a similar physics to a large extent.

Introduction

In an isotropic system with a single locally conserved quantity, the energy, Fourier's law asserts that the heat flux J is proportional to the temperature gradient:

J = κ(β)∇β (1)
where β = (k B T ) -1 denotes the inverse temperature and where the scalar quantity κ(β) is the thermal conductivity. 1 This is a macroscopic law that becomes valid at the scales where the system can locally be described by a Gibbs state at inverse temperature β. Based on general thermodynamical arguments, one concludes that κ(β) ≥ 0, since the entropy production rate γ s is given by

γ s = ∇β • J = ∇β • κ(β)∇β ≥ 0.
This is however the only general property of κ(β) that can be derived from entropic considerations. As clearly realized by Green [?] and Kubo [?, ?], the value of the conductivity does depend on the microscopic dynamics through some dynamical correlations. It may be explicitly expressed as

κ(β) = 2 ∞ 0 dt dx J(x, t) • J(0, 0) β (2)
where • β denotes the equilibrium state at inverse temperature β, and where J represents now the current on the microscopic scale. The current-current correlator featuring in (2) may display a wealth of interesting phenomena that cannot be diagnosed by inspection of the equilibrium states alone.

In this review, we consider systems where diffusive transport on long space-time scales is expected, and we will mostly take this for granted. Instead, we focus on intermediate scales and we analyze the finite-time effects that determine the value of the conductivity, as predicted by (2). The systems that we will analyze are chains of classical or quantum oscillators with a Hamiltonian of the form H = H 0 + gH 1 . At g = 0, all oscillators are uncoupled and no conduction is possible. Moreover, we assume that the dynamics of each individual oscillator is integrable and anharmonic. For a classical oscillator, anharmonicity means that the frequency genuinly varies with energy, and for a quantum oscillator, it means that the level spacing is non-constant, see below for concrete examples. When g > 0, the Hamiltonian H 1 provides a coupling between the oscillators, but the anharmonicity prevents efficient transport of energy due to frequency mismatch, except in rare places called resonances where there is accidentally no mismatch. These resonances play a crucial role in the energy transfer, since they may host chaotic spots that will eventually break the quasi-periodic behavior and ensure the validity of Fourier's law (1) in the long run.

The analysis discussed below yields a completely different picture for transport than the one involving weakly interacting phonons, that prevails near the harmonic (ballistic) limit [?, ?, ?]. Indeed, close to g = 0, our systems exhibit several characteristics of dynamical glasses: strongly intermittent dynamics with only a very low density of active spots, constrained motion of the chaotic spots due to the random energy landscape, and drastic fall-off of their mobility as one approaches the limit g = 0 (which translates in most systems of interest into some limit for the temperature or some other thermodynamical parameter). This behavior was unveiled in [?, ?, ?, ?, ?, ?, ?, ?, ?] and was recently checked through a detailed numerical analysis in [?, ?]. Hence, the poor transport properties of systems close to the atomic limit is reasonably understood. Yet, the theory still lacks testable quantitative predictions for the time scales on which transport eventually sets in, despite the detailed analysis in [?] for the disordered discrete non-linear Schrödinger (DDNLS) chain, and the (related) one in [?] for the translation invariant rotor chain.

That frequency mismatch inhibits the transfer of energy is by itself a well-known physical phenomenon. In classical few-body Hamiltonian systems, the KAM theorem shows that, at small coupling, a large proportion of the phase space is occupied by invariant tori where the motion is quasi-periodic, see e.g. [?], while Nekhoroshev's estimates guarantee the conservation up to stretched exponentially long times in the perturbative parameter of the uncoupled actions, see e.g. [?]. In quantum mechanics, disorder may strongly suppress quantum tunneling and lead to the localization in space of the wave packet, as first discovered by Anderson [?]. These phenomena are rather fragile though, and they are not generally expected to survive in the presence of dissipative effects, e.g. coupling to an environment. Therefore, if we consider extended systems at positive values of the thermodynamical parameters, it is a priori all but clear whether and how those phenomena can persist. A careful investigation started in the mid 2000s for disordered quantum systems, opening the field of many-body localization (MBL) [?, ?, ?]. Remarkably, the localization in some systems was found to be so robust that it persists in the thermodynamic limit [?, ?, ?]. Such systems remain however rather scarce and, as we stressed above, we will deal in this review with a much broader class of systems where any form of integrability is eventually washed out by dissipative effects.

The rest of this review is organized as follows. For convenience, we restrict ourselves to onedimensional systems, but none of the considerations developed here depends on the dimension. In particular, poor transport is not due to the possible presence of bottlenecks in the chain [?, ?, ?], though these effects may be present and constitute an additional slowing down factor. The three next sections deal with classical systems. In Section 2, we study the rotor chain, or classical Josephson junction chain, as a generic example. Generalizations, together with some related questions, are discussed in Section 3. In Section 4, we make the connection with quenched disordered systems. We discuss quantum systems in Section 5.

The rotor chain

Let us consider the classical Hamiltonian

H(θ, ω) = m 2 L x=1 ω 2 x + g L x=1 (1 -cos(θ x+1 -θ x )) (3)
where L is the length of the system, where θ = (θ x ) 1≤x≤L with θ x ∈ R/2πZ are angles and where ω = (ω x ) 1≤x≤L with ω x ∈ R are angular velocities. For convenience, we assume periodic boundary conditions θ L+1 = θ 1 , and we take units of mass such that m = 1. This system has two locally conserved quantities: the energy H and the angular velocity Ω = L x=1 ω x . Let us consider the system at equilibrium at inverse temperature β and zero average velocity, as can always be assumed without loss of generality by performing a Galilean boost. The variables θ and ω are conjugated and the dynamics is given by Hamilton's equations:

θx = ω x , ωx = g (sin(θ x+1 -θ x ) -sin(θ x -θ x-1 )) (4)
In the limit g = 0, the system reduces to independent oscillators that are trivially integrable: the angular velocity of each of them is conserved, ω x (t) = ω x (0), and the angles evolve as θ x (t) = θ x (0) + ω x (0)t. Let us next consider g > 0 small enough such that βg 1. Though the angular velocities now evolve with time, their statistics remains unchanged in the Gibbs state: ω x are independent and identically distributed (iid) Gaussian random variables with mean 0 and variance β -1 . Moreover, in first order approximation in g, the evolution of ω x is given by

ω x (t) ω x (0) + g t 0 ds (sin(∆ω x+1 (0)s) -sin(∆ω x (0)s)) (5)
with ∆ω x = ω x -ω x-1 . The integral featuring in ( 5) is oscillatory and averages to 0, so that it yields at all times a correction for the angles with typical value of order βg. It is however crucial to realize the possibility of atypical behavior: First order resonances occur when one of the differences ∆ω x+1 (0) or ∆ω x (0) is atypically small, i.e. of order g 1/2 or smaller. In this specific system, due to the presence of a second conserved quantity, two coupled oscillators do still form together an integrable system (that can be mapped to a pendulum) and the above resonances cannot by themselves constitute a mechanism to transport the energy. However, the dynamics of three consecutive oscillators with their angular velocity in a band of width g 1/2 is chaotic. Such a configuration is called a chaotic triple, and it occurs with a probability βg. See [?] for a deeper discussion on the origin of chaos in Hamiltonian systems, and [?, ?] for specific applications relevant here.

From this crude first order analysis, one concludes already that the dynamics displays an intermittent behavior in the regime βg 1 (at least as long as first order corrections dominate, but we will see that the conclusion extends to much longer times): Most oscillators display a quasi-periodic behavior while rare resonant triples act as heat baths and redistribute the energy into the chain. Intermittency is illustrated on Fig. 1, where three snapshots of the dynamics at different times are shown. 2 See also Fig. 4 in [?].

Let us next analyze how the chaotic triples may induce transport. Since these spots occur with a probability βg, a contribution of order (βg) 2 to the conductivity could be a first naive guess, but it turns out to be much smaller. To see it, following [?], we contemplate two mechanisms for the transfer of energy and, incidentally, we notice already that each of them will play a distinguished role in different quantum systems, see Section 5. First, fixed spots may exchange energy among themselves. However, since they are far appart and since the material in between is assumed to contain no chaotic region, this will require higher order processes, involving typically (βg) -1 oscillators. Second, spots may redistribute the energy in the oscillators in their vicinity and start moving across the chain, through the formation and the destruction of resonances. However, the process of transferring energy between the spot and a site at a relative angular velocity ∆ω compared to the average angular velocity of the spot is reduced by a factor e -∆ω/g 1/2 ∼ e -(βg) -1/2 typically. This reduction also "explains" the validity of Nekhoroshev's estimates.

We came thus to the conclusion that there is no transport due to first order effects in βg (including first order resonant spots). It is possible to undertake a systematic order-by-order analysis [?, ?] of higher order analogues of (5). The conclusion remains always the same: there is no transport at any given order, but mechanisms are identified that will lead eventually to the transfer of energy on longer time scales. By an optimization procedure first carried over in [?] for the DDNLS chain, one comes to a prediction for the rate γ at which energy is transferred:

γ ∼ g 1/2 e -c(log(βg)) 2 , c > 0 , (6)
see also [?]. As it turns out, the most efficient channel for transport is the second one described above, i.e. chaotic spots move around with a rate proportional to γ given by ( 6) and they modify the energy landscape in their vicinity. This mechanism, as well as the very small value of γ, properly reveal the glassy dynamics of the rotor chain in the regime βg 1. It would be an overkill to take the analytical expression (6) literally, since its derivation involves several approximations and since it would be exceedingly hard to validate it numerically, let alone experimentally. It has nevertheless the virtue to highlighting two robust pieces of information. First, at finite volume, the rate γ would be a stretched exponential in the inverse coupling parameter, as follows from Nekhoroshev's estimates [?]. Rigorous upper bounds at finite volume yield a power in the stretched exponential decaying as 1/L as L → ∞, and the expression (6) confirms that this power indeed vanishes in the thermodynamic limit. We notice in passing the rigorous result in [?]: An upper bound for the rate γ for an initially localized wave packet in an infinite lattice (i.e. at zero density) is derived and it basically coincides with (6).

Second, (6) shows that the rate γ decays faster than any polynomial in βg in the limit βg → 0. This leads to the claim that, in this limit, the conductivity is non-perturbative or non-analytic at origin as a function of βg, see [?, ?]. This fact is supported by several mathematical results, discussed below, and should be kept in mind when interpreting numerical results. As an example, let us report on numerical simulations of a simple quantity aimed to capture the relaxation time in the system:

ϕ(g, β, t) = lim L→∞ β 2L L x=1 (ω x (t) -ω x (0)) 2 β .
Clearly, ϕ(g, β, 0) = 0 and we expect that ϕ(g, β, t) → 1 as t → ∞. Thanks to the scaling relation ϕ(g, β, t) = ϕ(g/α 2 , βα 2 , αt) valid for any α > 0,3 it suffices to analyze ϕ(g, β, t) at a given value of β and vary g (and we further drop β from the notations). In Fig. 2(a), we show the the time evolution of ϕ for different values of g as a function of the rescaled time τ (g) = t/t 0 (g) with t 0 (g) = min{t : ϕ(g, t) = 0.9}. Under this rescaling, we observe that the data fit very well with each others, despite some very small but neat deviation for the smallest value g = 0.1, and we conclude that t 0 (g) can serve as a meaningful indicator of the relaxation time of the dynamics. 4The dependence of t 0 on g is shown on Fig. 2(b) and a fit with a power law yields t 0 (g) ∼ g -5 . 5 This is in line with recent numerical results in [?] where a similar parameter measuring the ergodization time is found to scale as g -6.5 . One can interpret these observed high powers as a confirmation of the drastic increase of time scales in (6) as g → 0. If taken literally, these power laws would actually contradict ( 6) and the non-perturbative nature of transport. However, since the bounds on the perturbative contributions to transport are actually mathematically rigorous (see below), we think that what is at play here is the numerical difficulty of distinguishing a high power from a superpolynomial function. Stated differently, it is hard to know whether the proper asymptotic regime is reached. Finally, we would like to stress that the non-analytical behavior of the conductivity as βg approaches 0 can be justified mathematically to a large extent, though not entirely. To get a tractable mathematical problem, it is convenient to introduce some noise in the system, following a strategy developed in [?] in the context of heat conduction. The presence of noise is unavoidable in any real or even numerical experiment and, since the absence of transport is usually associated to integrable effects, we expect the noise to enhance the dissipation of the energy. The upper bound below on the conductivity for the noisy system is thus expected to hold for the isolated system too.

We introduce a noise that conserves the total energy (but incidentally not the angular velocity), as it simplifies the theoretical analysis of the conductivity. Practically, the noise flips independently the velocity of the oscillators at a rate λ. Formally, the full dynamics is now described by the stochastic generator L = A + λS, where

A = ω • ∇ θ -∇ θ H • ∇ ω is the usual Liouville operator and Sf (θ, ω) = L x=1 f (θ, ω 1 , . . . , -ω x , . . . , ω L ) -f (θ, ω) . (7)
Let us assume that βg 1 as before, and let us set

λ = g 1/2 (βg) r for r 1 (8)
i.e. λ is a high power of βg. The mathematical result shown in [?] claims that, for fixed β > 0 and for any r > 0 there exists a constant C r (β) < +∞ so that

κ(g, β, λ) ≤ C r (β)g 2 λ. (9)
The interpretation of this bound is clear: There is no mechanism for transport that is stronger than the added noise with strength λ. Back to the system at λ = 0. It is convenient to deal with the reduced dimensionless conductivity κ(g, β) = g -1/2 β 2 κ(g, β) since scaling relations used above yield here that κ depends only on g, β through the product βg. Loosing strict mathematical rigor, we infer from (9) that κ(βg) ≤ C r (βg) r for all r > 0.

(10)

The only way that the latter bound would fail is if there was actually a transport mechanism that was destroyed by the added noise of strength λ, something that we find very hard to imagine in chaotic diffusive systems6 .

Translation invariant classical systems

The phenomenology discussed above for the case of a chain of rotors applies, or is expected to do, to a large class of translation invariant classical chains of oscillators. In [?], the discrete non-linear Schrödinger (DNLS) chain

H(ψ, ψ) = U L x=1 |ψ x | 4 + g L x=1 |ψ x+1 -ψ x | 2 , (11)
with the dynamics i ψ = ∇ ψ H, is considered in the regime g/U 1 or analogously for large value of the action density α = L x=1 |ψ x | 2 /L. Results comparable to ( 9) and (10) are found. Since both the rotor chain (3) and the DNLS chain (11) do have two conserved quantities, it is worth stressing that it did not play any role in the above discussion. In particular, the glassy behavior that we discuss here should be distinguished from the one that may arise due to the frustration of two conserved quantities, even though they have a common dynamical origin. In the DNLS chain (11), if the energy density ε = H/L exceeds a certain threshold for a given value of the action density α, the system enters in an effective negative temperature phase and forms a macroscopic condensate, or dynamical breather [?, ?, ?, ?, ?]. A breather turns out to be stable on extremely long time scales, due to a very strong frequency mismatch, and leads to a very slow approach towards equilibrium. In contrast to the phenomenology considered in this review, the formation of such breathers in the effective negative temperature phase is driven by entropic forces. It does not need to entail slow evolution for all observables of interest, and in particular the thermal conductivity in this phase is expected to behave just as the conductivity in the usual infinite temperature phase (i.e. a Gibbs state with no constraint on the energy), at least in d > 1 where no bottleneck effects are present.

More generally, we may consider oscillator chains of the form

H(p, q) = L x=1 p 2 x 2m + U (q x ) + V (q x+1 -q x ) ( 12
)
where the momenta p and the positions q are conjugated variables, and e.g. q L+1 q L (free boundary conditions). We expect bad transport properties at high temperature if the on-site potential U grows faster at infinity than the coupling potential V :

U (q) ∼ q a and V (q) ∼ q b as q → ∞ with a > b ≥ 2 (13)
where the condition a > 2 guarantees that the on-site frequency grows with the energy. To our knowledge, transport in the chain ( 12) under the conditions ( 13) has not been properly investigated so far, but several indirect studies are worth pointing out.

One way to understand local dissipative properties of extended systems is to take a small chain and connect it at the boundary to perfect heat baths, modeled e.g. by Langevin stochastic thermostats. In [?], the chain ( 12) is considered for polynomial potentials with a > b = 2 in ( 13). An effective model for the dissipation of the energy into the baths is derived, assuming that there is a single dynamical breather, i.e. a single oscillator at atypically high energy, initially in the middle of the chain:

dH dt ∼ -H -f (L) with f (L) = 1 - 2 a (L -1) -1 . ( 14
)
Numerical tests show a very good agreement between the prediction from ( 14) and the true behavior of the system. It should be noticed that, despite the appearance of the length L in this expression, it is only valid in the regime where L is fixed first and the energy of the breather is then taken to be very large. A similar analysis is performed for the rotor chain (3) in [?] coupled on one side to a thermostat at zero temperature, inducing only friction, so that all rotors eventually come to rest. A lower bound for the dissipation rate, of a similar kind as ( 14), is derived mathematically and shown numerically to be optimal.

If the chain ( 12) is connected at the two ends to reservoirs at different temperatures, it is expected to relax to a non-equilibrium stationary state (NESS). Moreover, for a fixed length, one would typically expect the relaxation to be exponentially fast. This is indeed what happens if b < a, i.e. if the condition ( 13) is violated so that the pinning becomes irrelevant at high energies, see e.g. [?] and references therein. Instead, when (13) holds and when the convergence to the NESS can be established, it is shown in [?] that the decay towards the NESS only follows a stretched exponential, due to the presence of breathers in the system. The same conclusion holds for chains of three and four rotors [?, ?, ?] (the case of more than four rotors has not been investigated mathematically but a stretched exponential convergence is expected there as well).

Quenched disordered systems

In the rotor chain discussed in Section 2, the frequency mismatch between nearby rotors was due to thermal fluctuations. Another way to realize the same physical effect is to introduce quenched disorder, as e.g. in the classical Hamiltonian

H(p, q) = L x=1 p 2 x 2m + mω 2 x q 2 x 2 + g 0 (q x -q x+1 ) 2 + g q 4 x 4 (15)
with free boundary conditions (q L+1 = q L ) and where the frequencies ω x of the on-site harmonic oscillators are taken to be strictly positive iid random variables with standard deviation W > 0.

For this system, it may seem that the frequency mismatch is much more robust than in the rotor chain (3) since the frequencies ω x do not evolve with time. This is indeed so in the purely harmonic case g = 0, where the Hamiltonian (15) turns out to be an Anderson insulator in disguise. In particular, in d = 1, all eigenmodes of the chain are localized for any value of g 0 , implying a vanishing thermal conductivity [?] and a strong breakdown of local equilibrium [?], see also [?] and references theirein. However, the situation changes in the presence of anharmonic interactions, i.e. g > 0, that provide effective couplings among the localized modes. It was observed numerically in [?] that "even a small amount of anharmonicity leads to a J ∼ 1/L dependence, implying diffusive transport of energy", see also [?, ?]. Normal transport at any non-zero value of the coupling g was confirmed in [?] for a classical spin chain, where moreover a drastic reduction of the conductivity as g → 0 is observed, "consistent with the possibility that the transport is actually essentially non-perturbative in g". This was justified mathematically in [?] in very much the same way as for the rotor chain discussed in Section 2: Fixing the disorder strength W and taking here g/βg 2 0 instead of gβ as the perturbative parameter, the noise ( 7) is added to the dynamics with a very slow rate (8), leading to expressions analogous to (9) and (10) for the conductivity (see [?] for scaling relations of the thermal conductivity with the various parameters of the model in the chain (15)).

However, the origin of transport in the chain ( 15) is harder to elucidate than in the rotor chain. Indeed, there are now two types of chaotic spots: spots induced by fluctuations of the quenched disorder or by thermal fluctuations as in the rotor chain. The first one are abundant but they are obviously unable to travel through the chain. The second one can travel but require much more atypical thermal fluctuations than in the rotor chain to be formed, due to the presence of random frequencies. The detailed analysis carried in [?] for the related DDNLS chain shows that the rate of dissipation is still of a strength comparable to (6), i.e. smaller than any power law, but larger than any exponential. More important, it turns out that here as well, thermal fluctuations furnish the fastest channel for transport. If this theory is correct, we are thus led to the conclusion that the presence of quenched disorder is eventually irrelevant and that the dynamical behavior becomes glassy in the chain (15) as well, in very much the same sense as for the translation invariant systems.

Finally, we would like to mention some puzzling observations on the spreading of an initially localized wave packet in a disordered medium where all other oscillators are initially at rest, as they might hint at the fact that some important aspect of the problem has been overlooked in the discussion so far. As is it the case for the chain in (15), we consider systems where the limit g → 0 may be equivalently realized by taking the limit β → ∞, which is relevant here since, if any local equilibrium sets-in, the local temperature must go to 0. This set-up has been studied intensively. Theoretical [?, ?, ?] and mathematical [?, ?, ?, ?, ?] studies all yield results on the slow, essentially non-perturbative, spreading of the wave packet as a function of time, while the evolution is fit with a power law in most numerical studies, see e.g. [?, ?, ?, ?, ?, ?, ?]. Mathematical results are however hard to compare formally with numerical observations, since the former deal either with a transient regime or impose strong constraints on the parameters of the model. To see that these observations are conflicting with the arguments developed here, let us focus on the results in [?, ?], where the second moment m 2 of the wave packet is found to scale as t α with α = 1/3. Further, it is argued in [?, ?] that this exponent describes the proper asymptotic long-time regime, that the value α = 1/3 is universal, and that the "spreading is due to corresponding weak chaos inside the packet, which slowly heats the cold exterior" [?]. If the packet indeed spreads through local thermalization, then the observed scaling of m 2 is incompatible with the fact that the conductivity decays faster than any power law in the temperature.

Quantum systems

The main reason why quantum systems may behave differently from classical ones is that a few degrees of freedom cannot be properly chaotic, i.e. the spectrum of a finite system remains discrete and cannot behave as a perfect bath. And indeed, as pointed out in the introduction, a specific class of systems was found to host an MBL phase, characterized by the total absence of transport of the conserved quantities and an emergent integrability, see e.g. the reviews [?, ?, ?]. Systems known to display genuine MBL are one-dimensional quenched disordered lattice Hamiltonians with a finite Hilbert space per site, such as chains of spins or fermions on a lattice. Extensions have naturally been investigated and in particular, the possibility of MBL in clean, translation invariant systems is of interest. Investigated cases include: systems with long range anharmonic interactions [?] (even much before MBL was coined), long lived metastable states for lattice bosons [?], two species models [?, ?, ?], Bose-Hubbard Hamiltonian at high density of particles [?, ?], multi-lane models [?], quantized classical glasses [?], and spin chains with strongly anharmonic interactions [?]. As it turns out, these have only be found to exhibit slow dynamics but no true MBL phase [?, ?, ?].

To get one concrete example, let us consider the Bose-Hubbard Hamiltonian

H = U L x=1 n 2 x + g L-1 x=1 (a † x a x+1 + a † x+1 a x ) (16)
with n x = a † x a x , that can be regarded as the exact quantum analog of the DNLS chain in (11). Results on the non-perturbative character of the conductivity in the limit where the density of particles ρ → ∞ have been derived in [?], improving on the techniques first developed in [?]. For the quantum chain (16), there is no obvious way to introduce noise in the system and instead, to obtain a mathematical statement, the time-integral defining the conductivity (2) is truncated to a high power in ρ. The mathematical result in [?] supports the claim that

ρ r κ(ρ) → 0 as ρ → ∞ for all r ≥ 0,
for fixed values of the temperature and the parameters of the model g, U .

Let us next analyze why this system is not expected to be genuinely localized (in particular why the conductivity does not vanish exactly starting from some density ρ). Obviously, one may first argue that this is just impossible since the Hamiltonian in (16) approaches its classical limit (11) as ρ → ∞. This is true but this difficulty is superficial and can be circumvented. For example, in the spirit of [?], one may impose a cut-off on the total number of bosons per site, and consider the system in the regime g/U 1 rather than in the limit ρ 1. A similar set-up is considered in [?, ?, ?].

This said, following [?, ?], let us now come to a reason for delocalization that applies to all models quoted above, but let us stick to the Hamiltonian (16) for commodity. Due to thermal fluctuations, one encounters regions where the density of particles is anomalously low, i.e. much lower than ρ. There, the system is expected to be thermal and these spots act thus as (imperfect) baths. Even though they are imperfect, it should be realized that the level spacing decreases exponentially as the size of the rare regions grows, hence a large system contains unavoidably spots that can be considered effectively as good baths, in the sense that the sites in their vicinity can be thermalized. Therefore, by heating up near spins, the spot can migrate into the chain and modify the energy landscape along the way. We come thus to the conclusion that the main mechanism for energy transfer in strongly anharmonic quantum chains near the integrable limit is properly the same as for classical systems, even though more space than two or three sites may be needed to realize a good enough effective bath.

Finally, it is interesting to notice that rare regions may also have deep destabilizing effects on quenched disorder spin chains. First, the possibility of hot mobile bubbles, as described above, was shown to imply the absence of many-body mobility edges [?], i.e. that the presence of thermal states at a given value of the thermodynamical parameters implies that states are thermal at all values of these parameters. Second, while for classical systems we concluded in Section 4 that the presence of chaotic spots due to fluctuations of the quenched disorder did not play a major role in the transport of the energy, their presence become relevant for quantum systems: They destabilize the MBL phase in d > 1 [?, ?] and are expected to drive the system to the thermal phase if the bare localization length exceeds some threshold value in d = 1 [?, ?].

Figure 1 :

 1 Figure 1: Three snapshots of the evolution of the angular velocities ω x for x = 1 to x = 40 for the Hamiltonian (3) at thermal equilibrium, with L = 64, m = 1, β = 1 and g = 0.07. Vertical tick marks below correspond to the values ω x = 0 for each oscillator.

Figure 2 :

 2 Figure 2: Panel (a). Time evolution of ϕ as a function of the rescaled time τ (g) for various values of g (see main text), for L = 64, β = 1 and m = 1. Panel (b). Relaxation time t 0 (g) for the same system (see main text). The dotted line with slope -5 serves as a guide for the eye. Average over 8000 initial configurations at least for each value of g.

It is more conventional, but slightly less convenient here, to define the thermal conductivity through the relation J = -κ(T )∇T with k B κ(T ) = β

κ(β).

In our simulations, the system is first put in thermal equilibrium by acting on each coupled oscillators with a Langevin heat bath at inverse temperature β over a time t th = 1000, and the isolated dynamics is then simulated with a Störmer-Verlet scheme with an elementary time step ∆t = 0.001 for Fig.1and ∆t = 0.01 for Fig.2.

This follows from the scaling relations θ(t, g; θ o , ω o ) = θ(αt, g/α 2 ; θ o , ω o /α) and ω(t, g; θ o , ω o ) = αω(αt, g/α 2 ; θ o , ω o /α), with (θ o , ω o ) the initial conditions and α > 0, that are derived from (4).

The data strongly suggest that ϕ grows monotonically as a function of t. To determine t 0 , we have averaged the values of ϕ in small intervals of time so as to get smooth monotonic curves.

We did not find any dependence on the system size L and this was also not expected: Since the angular velocity is conserved, we expect that 1 -ϕ(g, t) ∼ c(g)t -1/2 as L → ∞.

This would, of course, make sense in integrable ballistic systems, but that is not the focus here.

Acknowledgements. This work of F.H. was partially supported by the grants ANR-15-CE40-0020-01 LSD and ANR-14-CE25-0011 EDNHS of the French National Research Agency (ANR). W.D.R. acknowledges the support of the Flemish Research Fund FWO under grant G076216N.


Daniel F Araújo

Emmanuel Ponzevera

Nicolas Briant

Joël Knoery

Sandrine Bruzac

Teddy Sireau

Christophe Brach-Papa

Copper, zinc and lead isotope signatures of sediments from a mediterranean coastal bay impacted by naval activities and urban sources

Keywords: metal isotopes, isotope tracers, geochemical isotopes, environmental pollution, coastal pollution, Toulon Bay

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The trace metal contamination of coastal systems, notably by copper (Cu), zinc (Zn) and lead (Pb), is a global environmental problem that imposes major risks for human health and ecological damages [START_REF] Barletta | Distribution, sources and consequences of nutrients, persistent organic pollutants, metals and microplastics in South American estuaries[END_REF][START_REF] De Souza Machado | Metal fate and effects in estuaries: A review and conceptual model for better understanding of toxicity[END_REF][START_REF] Naser | Assessment and management of heavy metal pollution in the marine environment of the Arabian Gulf: A review[END_REF].

Toulon bay (NW Mediterranean Sea) is an emblematic coastal environment severely impacted by a multi-metallic contamination for Cu, Zn, Pb, As, and Hg originating from legacy to present recent naval and seashore human activities (marinas, industry, ship-building and maintenance, Navy). While a large load of these contaminants was released into the bay during the Second World War (WWII) as result of the French Navy fleet bombardment and scuttling [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF]Grasset, 2011), it has been reported trends of increasing on metal releasing by diffuse sources, including urban effluents, industry and anti-fouling paints [START_REF] Coclet | Trace metal contamination as a toxic and structuring factor impacting ultraphytoplankton communities in a multicontaminated Mediterranean coastal area[END_REF][START_REF] Cossa | A Michaelis-Menten type equation for describing methylmercury dependence on inorganic mercury in aquatic sediments[END_REF]Dang et al., 2015a,b;[START_REF] Jean | Protein expression from zooplankton communities in a metal contaminated NW mediterranean coastal ecosystem[END_REF][START_REF] Pougnet | Sources and historical record of tin and butyl-tin species in a Mediterranean bay (Toulon Bay, France)[END_REF][START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF].

Metal inputs of Cu, Zn and Pb into the bay since-WWII period are estimated to be about 595, 1091 and 817 tons respectively [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF]. Most part of this stock (about 80%) rests confined in sediments along the military and civils ports of the bay, favored by the long water residence time and low hydrodynamic energy [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF]. However, the rest of this metal inventory (~ 20%) is accessible to resuspension and transport toward the open sea associated to sediment particles (Dang et al., 2015a). This availability imposes major ecological and health risks to marine trophic webs and to local aquaculture activities, like as fishing and aquaculture farming (Dang et al., 2015b;[START_REF] Rossi | In situ heavy metals (copper, lead and cadmium) in different plankton compartments and suspended particulate matter in two coupled Mediterranean coastal ecosystems (Toulon Bay, France)[END_REF]. Thus, understanding the origin of metal sources, pathways and sinks in this ecosystem are of upmost priority to develop effective environmental management, policy strategies, and monitoring programs.

To get a better appreciation of spatial distribution of sources and dynamics of metal contaminants, we combined complementary information from the stable isotopes of copper ( 63 Cu, 65 Cu), zinc ( 64 Zn, 66 Zn, 67 Zn, 68 Zn) and radiogenic Pb isotopes ( 206 Pb, 207 Pb, 208 Pb) in conjunction with a high-resolution surface sediment sampling in Toulon bay. The variability of Pb isotope compositions in environmental samples is induced by the radioactive decay of uranium (U) and thorium (Th) radionuclides and depends on the initial proportion between the parental and radiogenic nuclides [START_REF] Komárek | Lead isotopes in environmental sciences: A review[END_REF]. The magnitude of lead isotope fractionation due to natural and anthropic processes is negligible in comparison to variations induced by U and Th decay. Thus, Pb isotope ratios (eg., 206 Pb/ 207 Pb) are conservative through biogeochemical processes [START_REF] Komárek | Lead isotopes in environmental sciences: A review[END_REF]. Therefore, their ratios are convenient environmental tracers in the diverse or multiple pollution contexts like the urban atmosphere [START_REF] Souto-Oliveira | Multi-isotope approach of Pb, Cu and Zn in urban aerosols and anthropogenic sources improves tracing of the atmospheric pollutant sources in megacities[END_REF], 2018), mining operations and metallurgy [START_REF] Dolgopolova | Use of isotope ratios to assess sources of Pb and Zn dispersed in the environment during mining and ore processing within the Orlovka-Spokoinoe mining site (Russia)[END_REF], urban sewage, and agricultural practices [START_REF] Charalampides | Sr and Pb isotopes as environmental indicators in environmental studies[END_REF]. Zinc isotopes fractionate significantly during redox and evaporation processes during ore refining (smelting and electroplating industry ), and coal combustion [START_REF] Gonzalez | Zinc Isotope Variability in Three Coal-Fired Power Plants: A Predictive Model for Determining Isotopic Fractionation during Combustion[END_REF][START_REF] Kavner | Redox-driven stable isotope fractionation in transition metals: Application to Zn electroplating[END_REF]Yin et al., 2015). This results in zinc manufactured products and byproducts being isotopically distinct from their natural compartments [START_REF] Gonzalez | Zinc Isotope Variability in Three Coal-Fired Power Plants: A Predictive Model for Determining Isotopic Fractionation during Combustion[END_REF][START_REF] Shiel | Evaluation of zinc, cadmium and lead isotope fractionation during smelting and refining[END_REF][START_REF] Sivry | Zn isotopes as tracers of anthropogenic pollution from Zn-ore smelters The Riou Mort-Lot River system[END_REF], allowing to trace Zn anthropogenic inputs and their transfers in the sediment-water-biota interfaces [START_REF] Mar Da Costa | Biogeochemical mechanisms controlling trophic state and micropollutant concentrations in a tropical artificial lake[END_REF]2017a,b;2018, 2019). Unlike Zn, Cu does not fractionate during smelting ore refining due its high boiling point (2562 °C). Manufactured, slag and other metallurgical by-products inherit the isotope signature of Cu ore concentrates (Šillerová et al., 2017). As Cu isotope compositions of minerals formed in high temperature (>250 °C)such as chalcopyrite, CuFeS 2 -are different from low temperature (<100°C) minerals -such as chalcocite, Cu 2 S-, Cu has become a promising tool to trace origins of archeological materials and archeometallurgical processes [START_REF] Jansen | The potential of stable Cu isotopes for the identification of Bronze Age ore mineral sources from Cyprus and Faynan: results from Uluburun and Khirbat Hamra Ifdan[END_REF]Klein et al., 2009). Since anthropogenic Cu may derive either from sulfides or oxide ore deposits, Cu isotope systematics can help discriminate the different anthropogenic and natural sources. In environmental studies, Cu isotopes have been a useful tracer in urban aerosol samples [START_REF] Souto-Oliveira | Multi-isotope approach of Pb, Cu and Zn in urban aerosols and anthropogenic sources improves tracing of the atmospheric pollutant sources in megacities[END_REF], 2018) and in surrounding soils impacted by metallurgical (Mihaljevič et al., 2019;[START_REF] Šillerová | Stable isotope tracing of Ni and Cu pollution in North-East Norway: Potentials and drawbacks[END_REF] and vineyard activities [START_REF] Babcsányi | Copper in soil fractions and runoff in a vineyard catchment: Insights from copper stable isotopes[END_REF][START_REF] El Azzi | Origin and fate of copper in a small Mediterranean vineyard catchment: New insights from combined chemical extraction and δ65Cu isotopic composition[END_REF], however their applicability in complex and reactive aquatic systems has been more challenging [START_REF] Araújo | Assessment of the metal contamination evolution in the Loire estuary using Cu and Zn stable isotopes and geochemical data in sediments[END_REF][START_REF] Petit | Anthropogenic sources and biogeochemical reactivity of particulate and dissolved Cu isotopes in the turbidity gradient of the Garonne River (France)[END_REF][START_REF] Viers | Are Cu isotopes a useful tool to trace metal sources and processes in acid mine drainage (AMD) context?[END_REF].

The particular feature of Toulon bay with anomalous metal enrichment in nearshore sediments contrasting against more pristine offshore sediments, provides a unique opportunity to test a multi-isotope approach (Cu-Zn-Pb) as an environmental forensic tool able to deconvolute spatial source distribution of metal and inferring dispersion mechanisms for these contaminants along the land-sea continuum.

Materials and methods

Study area and sampling

Toulon Bay is an urbanized French coastal ecosystem of approximately 600,000 inhabitants separated from the open sea to the East and South by the peninsula of Giens, and St.Mandrier, respectively (Tessier et al., 2011, Fig. 1). Two main urban rivers bring freshwater into this bay: the Las and Eygoutier rivers (Misson et al., 2016, Fig.1). A seawall built in the 19 th century nearly isolates two sub-basins labelled the "Small bay" (9.8 km 2 , semi-enclosed), on the northwestern side, and "Large bay" (42.2 km 2 ), on the southeastern part (Fig. 1). On the Small bay shore are located urbanized areas, an arsenal and military port, industrial, recreational and commercial ports and former shipyards (closed in 1989). The

Large bay is open to Mediterranean Sea to its east, and hosts a single military zone and several small beaches and tourist ports. A wastewater treatment plant releases treated effluent to the open sea close to sampling point 34 (Fig. 1), and a second in the western, outer part of the Bay (Cap Sicié, not visible in Fig. 1).

Fifty-five sampling sites of surface sediment (0-5 cm) were considered to this study (Fig. 1), accounting for twenty-four sites in the Small bay (Cs, Mis, Petite Mer, SMS and Stations 1-19) and twenty-eight in the Large bay (Stations 21-52). The samples were collected between November 2008 and June 2009 in the frame of the CARTOCHIM project driven by the PROTEE laboratory at the University of Toulon. They were characterized for their trace metal concentrations and Pb isotopes (Dang et al., 2015b;[START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF], Sn and butyl-Sn species in previous studies [START_REF] Pougnet | Sources and historical record of tin and butyl-tin species in a Mediterranean bay (Toulon Bay, France)[END_REF]. Here, new aliquots from the same sediment samples were dissolved for trace elements and Cu, Zn and Pb isotope analyses.

Sample preparation and instrumental analysis

All reagents, labware acid cleaning and solution dilutions for elemental and isotope analyses were performed using 18.2 MΩ.cm -1 H 2 O (Nanop System®) and ultra-pure acids (PlasmaPure Plus grade, SCP science®). Unsieved dry aliquots of sediments were digested in Teflon® bombs on a coated graphite block using multiple-step acid procedure with HF, HCl and HNO 3 . The final extract solution was split for subsequent elemental and isotope analyses.

Procedural blanks and international certified reference materials (MESS-3 and PACS-2 from NRC -CNR®) were processed in each batch of sediment samples and followed the complete procedure (digestion, elemental analysis, ion chromatography and isotope measurements).

Elemental concentrations of Cu, Zn, Pb, Cd, Ag, Ni, Co, Cr, V, Li, Mn, Al and Fe were determined by inductively coupled plasma mass spectrometry (Q-ICP-MS, iCAP Qc, Thermo Fisher Scientific) at the Laboratoire de Biogéochimie des Contaminants Métalliques (LBCM, Ifremer, Nantes). The accuracy, as average bias, was always within 10% of the certified values of both reference materials for considered elements. Lead isotope ratio ( 206 Pb/ 207 Pb and 208 Pb/ 206 Pb) measurements were performed with the same instrument. Mass bias and instrumental drift were corrected with a standard bracketing method using the NIST SRM-981 standard reference material. The 206 Pb/ 207 Pb and 208 Pb/ 206 Pb internal relative standard deviations average were of 0.20 ±0.07% and 0.47±0.17% (1s, n = 60), respectively, for unknown samples and sediment reference materials (MESS-3 and PACS-2). Lead isotope ratios obtained for the reference materials are included in the Table 1S (Supplementary Material).

For isotope analysis, Zn and Cu were isolated from matrix elements by ion-exchange chromatography technique adapted from protocol of [START_REF] Maréchal | Precise analysis of copper and zinc isotopic compositions by plasma-source mass spectrometry[END_REF]. The whole chromatography procedure was performed in an overpressured air-filtered cleanroom (Class 1000) at the Laboratoire de Biogéochimie des Contaminants Métalliques (LBCM, Ifremer, Nantes).The recovery yields (around 100 ±10%) were systematically checked by ICP-MS at the same laboratory. The total procedural blanks were below 1% of Zn and Cu mass content of the samples. Copper and Zinc isotopes were measured at the Pôle Spectrométrie Océan (Ifremer, Brest, France) using MC-ICP-MS (Neptune, Thermo Scientific). Samples were dissolved in 2% (v/v) HNO 3 and introduced at concentrations between 100 and 300 ng.g -1 with a stable introduction system (SIS: cyclonic spray chamber) and a low flow PFA nebulizer (50 µL.min -1 ). The final isotope compositions are expressed as δ 65 Cu and δ 66 Zn values in relation to NIST SRM 976 and to the "JMC-Lyon" (Johnson Matthey Company 3-0749-L), respectively, as follow:

( ( ) ( ) ) (Eq.1) ( ( ) ( ) ) (Eq.2)
Sediment reference materials measured three times along each analytical session were used for analytical quality control. The δ 65 Cu NIST and δ 66 Zn JMC average values obtained for marine sediments (MESS-3 and PACS-2) are reported in the Table 1S. External precision (2s) average obtained from two to three measurements of unknown samples and sediment certified materials was ±0.05‰ and ±0.06‰ for δ 65 Cu NIST and δ 66 Zn JMC, respectively. This corresponds to the uncertainty bars presented in the data plots (Fig. 5 and6).

Enrichment Factor

The Enrichment Factor (EF) of trace metals is an index used to estimate the degree of anthropogenic contamination. Dilution effects related to grain size and mineralogical heterogeneities are decreased by the normalization of a metal concentration in the sample to the concentration of a reference element (Chen et al., 2009). In this work, aluminum (Al) was used as reference element and EF values were calculated as follows:

( ) (Eq.3)
where "M" is the concentration of the metal of interest and "Al" is the reference element concentration, whereas "(M/Al) sample " represents the concentration ratio in the sample and "(M/Al) background " refers to the concentration ratio considered as free of anthropogenic influence, i.e., the representative of natural background. In this work, the sample station 31 was chosen as a representative natural sample based on its Pb isotope compositions ( 206 Pb/ 207 Pb = 1.1998, Table 1S), which is consistent with natural Pb isotope signatures referred in previous studies performed in Toulon bay ( 206 Pb/ 207 Pb natural = 1.195;Dang et al., 2015b). EF values were classed as suggested in other studies (Chen et al., 2007): EF < 1 indicates no enrichment; EF < 3 is minor enrichment; EF = 3-5 is moderate enrichment; EF = 5-10 is moderately severe enrichment; EF = 10-25 is severe enrichment; EF = 25-50 is very severe enrichment; and EF > 50 is extremely severe enrichment.

Statistics and geostatistical methods

All descriptive and multivariate statistics were performed using the SPSS® software v.18.0 (IBM®). Principal Component Analysis (PCA) was applied as exploratory tool to examine and identify the main interrelations among the isotope and elemental variables. PCA is useful for simplifying or extracting main dataset information by grouping inter-correlated variables in a smaller number of underlying factors, also known as principal components [START_REF] Meglen | Chemical multi-contamination drives benthic prokaryotic diversity in the anthropized Toulon Bay[END_REF]. Before run PCA, linearity of variables was checked by scatter plots, while sampling adequacy of the dataset and its suitableness for data reduction was validated with the Kaiser-Meyer-Olkin (KMO) index and Bartlett's test, respectively. PCA factor loads were calculated using the correlation matrix with varimax rotation. Factor loads express the influence of each variable in the principal components [START_REF] Field | Discovering Statistics Using IBM SPSS Statistics: North American Edition[END_REF]. The sample factor scores were calculated using the regression method to assess the influence of each component in individual samples [START_REF] Field | Discovering Statistics Using IBM SPSS Statistics: North American Edition[END_REF]Araújo et al., 2017c;). The principal components (PCs) that had eigenvalues higher than 1 were selected for discussion [START_REF] Mar Da Costa | Biogeochemical mechanisms controlling trophic state and micropollutant concentrations in a tropical artificial lake[END_REF][START_REF] Mulholland | Geological and anthropogenic influences on sediment metal composition in the upper Paracatu River Basin, Brazil[END_REF]. The communalities values of all variables equal or higher than 0.5 were used as thresholds for PCA results acceptance (Araújo et al., 2017c).

Interpolation maps were used for a better spatial visualization of elemental and isotope data and were performed using the Inverse Distance Weighting (IDW) method with ArcGIS® software v. 10.5(ESRI®). The IDW is the simplest interpolation method that uses the surrounding measured values to derive a prediction for an unmeasured location [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF].

Results and discussions

Enrichment factor values and spatial distribution of trace metals

All elemental isotope data for trace metals are presented in Table 1S (Supplementary Material). The calculated enrichment factors are shown in Table 2S (Supplementary Material). Copper, Zn, Pb, Ag and Cd present the highest enrichment factor values in sediments of Toulon bay (Fig. 2). In the Small bay, Zn, Ag, Pb reached severe enrichment levels (10 < EF < 25), while Cu and Cd presented a very severe enrichment (25 < EF < 50).

EF decrease significantly for these elements in the Large bay (Fig. 2), reaching levels of minor enrichment (EF < 3, for Zn, Ag, Pb, Cd) and moderate enrichment (3 < EF < 5, for Cu). It is noted in the Large bay two anomalous EF Cu values for the sampling stations 37 and 42 (Fig. 2).

The elements Li, V, Cr, Co and Ni display a homogenous distribution between the Small and Large bays, with low EF values (< 3). These low values indicate no enrichment or a minor enrichment, suggesting few or practically no anthropogenic influence (Fig. 2).

Overall, Mn also displayed low EF values (< 3) through the bay, but with a noticeable enrichment in the Large bay (average EF of 1.70) in relation to Small bay (average EF of 0.9).

This increase in Mn contents in surface sediments is apparent in the Fig. 2. The most exceptional sampling stations are 37, 47, 50, reaching the highest Mn EF values (3 < EF < 5, Fig. 2).

PCA analysis identified two main components, together explaining 82% of total variance of all dataset. The PCA component plots (Fig. 3a) provides a general pattern of metal to Mn(II), which is mobilized into the pore waters, and subsequently precipitating as Mnoxides or carbonates in the dominant oxic conditions of the Large bay (Araújo et al., 2017c;Dang et al., 2015a). This diagenetic Mn remobilization from the Small bay toward the Large bay explains the well-marked Mn enrichment zones between these systems (Fig. 2).

The plot of sample factor scores (Fig. 3b) allow verify the weight of each component obtained in the PCA analysis for individual sediment samples. For the surface sediments located in the Small bay, it is observed an association with elements of anthropogenic origin identified in the first component (Cu, Zn, Pb, Ag, Cd). On the other hand, both Small and Large bays present similar factor scores for the second component associated to natural elements (Li, Co, V, Cr and Ni). This reflects a homogenous natural metal baseline within the entire bay. Also, the associations and spatial distribution of trace metals match well with the spatial patterns of EF values observed in the maps of interpolation (Fig. 2). In summary, multivariate statistics and geostatistical approach shows a clear inter-bay contrast, highlighting the highest concentrations for Cu, Zn, Pb and Ag in the vicinity of civil and military ports with a strong punctual anomaly at the MIS station (collected inside the military zone). Moreover, an additional anthropogenic metal input coming from the basin drainage and delivered by the Las river water outlet is plausible.

In the following sections, Pb, Cu and Zn isotope systematics are discussed in order to characterize the natural and anthropogenic signatures and possibly to deconvolute potential multiple anthropogenic sources. Use of Cu and Zn stable isotopes as environmental tracers of surficial sediments requires that early diagenetic processes do not alter significantly source isotope records and that the sediment column remains undisturbed. Diagenetic reactions induced by microbial degradation of organic matter change redox conditions in sedimentary strata in the surrounding water, which in turn, promotes mineral dissolution or the formation of new mineral phases, including sulfides, carbonates and oxides [START_REF] Machado | Relation of Reactive Sulfides with Organic Carbon, Iron, and Manganese in Anaerobic Mangrove Sediments: Implications for Sediment Suitability to Trap Trace Metals[END_REF][START_REF] Mason | Trace Metals in Aquatic Systems: Mason/Trace Metals in Aquatic Systems[END_REF]. These mechanisms of adsorption of metals onto mineral surfaces, or in their inclusion in crystal lattices control their mobility and speciation and affect their isotope fractionation [START_REF] Dekov | Ferrihydrite precipitation in groundwater-fed river systems (Nete and Demer river basins, Belgium): Insights from a combined Fe-Zn-Sr-Nd-Pb-isotope study[END_REF][START_REF] Ducher | Equilibrium zinc isotope fractionation in Zn-bearing minerals from first-principles calculations[END_REF][START_REF] Jamieson-Hanes | Evaluating zinc isotope fractionation under sulfate reducing conditions using a flowthrough cell and in situ XAS analysis[END_REF][START_REF] Ratié | Nickel isotope fractionation during tropical weathering of ultramafic rocks[END_REF][START_REF] Veeramani | Zinc Isotope Fractionation as an Indicator of Geochemical Attenuation Processes[END_REF]. In Toulon bay, earlier examination of sediments and porewaters demonstrated increases of trace metals concentrations (Co, Cu, Pb) due to such processes mediated by Mn and Fe oxyhydroxides reduction under anoxic conditions (Dang et al., 2015a). However, the metal loads delivered from sediments to the overlying water column account for less than 1% of sediment pool, suggesting that changes of the Cu and Zn isotope compositions due to such biogeochemical processes are negligible. This is coherent with other studies that have also pointed out that Cu and Zn isotope compositions in bottom or particulate sediments are essentially controlled by source mixing processes (Araújo et al., 2017a;[START_REF] Babcsányi | Copper Stable Isotopes To Trace Copper Behavior in Wetland Systems[END_REF]Chen et al., 2009;Guinoiseau et al., 2018;[START_REF] Thapalia | Zinc Isotopic Signatures in Eight Lake Sediment Cores from Across the United States[END_REF][START_REF] Zhang | Historical records, distributions and sources of mercury and zinc in sediments of East China sea: Implication from stable isotopic compositions[END_REF]. The second criterion based on integrity of the stratigraphy is met since geochemical profiles of these trace metals in dated cores of Toulon bay were consistent with historical events of metal contamination indicating a good preservation of geochemical records [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF].

From an ecological perspective, understanding the diagenesis effects on elemental speciation and bioavailability at sediment-water interface are of prime importance to determine the bioaccumulation routes of metals in the trophic chain. This can occur directly from the sediment, as in the case of benthic organisms living in the sediments, or by the dissolved fraction of interstitial surface waters, in the case of primary producers (phytoplankton, macrophytes and bacteria) (Ahlf et al., 2009;[START_REF] Briant | Spatial and temporal distribution of mercury and methylmercury in bivalves from the French coastline[END_REF][START_REF] Chouvelon | Oligotrophy as a major driver of mercury bioaccumulation in medium-to high-trophic level consumers: A marine ecosystem-comparative study[END_REF][START_REF] Zhou | Biomonitoring: An appealing tool for assessment of metal pollution in the aquatic ecosystem[END_REF][START_REF] Zuykov | Bivalve mollusks in metal pollution studies: From bioaccumulation to biomonitoring[END_REF]. Then, a better understanding of the diagenesis influence on Cu and Zn isotope fractionation at sediment-water interface will be helpful to exploit the isotope records in biological tissues.

Pb isotope systematics

Bulk lead concentrations variated from 13 to 470 µg.g -1 , matching well with the previous concentration range reported by [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF]. The Small bay has an average Pb concentration of 144 µg.g -1 , while the Large bay has an almost threefold lower Pb concentration at around 49 µg.g -1 (Table 1).

Lead isotope ratios ( 206 Pb/ 207 Pb) are ranging from 1.1591 to 1.1998 in the entire sediment dataset. Averages ratios of 1.1664 ±0.0043 and 1.1763 ±0.0079 are observed in the Small and Large bays, respectively (Table 1). The lowest 206 Pb/ 207 Pb ratio values found in sediments from the Small bay are close to the anthropogenic signature value proposed by Dang et al., 2015 ( 206 Pb/ 207 Pb anthrop. ~ 1.166 ±0.003, Fig. 4), while sediments from the Large bay tend to approximate a more natural baseline ( 206 Pb/ 207 Pb natural ~ 1.18 -1.20, Fig. 4).

Plotting 206 Pb/ 207 Pb ratios vs. EF Pb values denotes a binary source mixing processes between natural and anthropogenic end members (Fig. 5a). The latter dominates in the Small bay and is attributed to multi-uses of Pb in the Navy and industrial activities (Dang et al., 2015b). The spatial interpolation map of 206 Pb/ 207 Pb ratios indicates a gradual attenuation of Pb contamination from the Small bay toward the open sea, as a result of transport and dilution by natural particles (Fig. 4).

Cu isotope systematics

Copper concentrations range from 4.4 to 666 µg. g -1 , while δ 65 Cu NIST ranged between -0.79 and +0.34‰ (Table 1). In the Small bay, the average Cu concentration reach 147 µg. g -1 and an average δ 65 Cu NIST of -0.17 ±0.19‰ (1s), while in the Large bay, Cu concentration average is lower (ca. 32 µg. g -1 ), and δ 65 Cu NIST tends to more positive and heavier values (+0.08 ±0.23‰, 1s).

The contrast between Cu isotopes fingerprints in sediments of Small and Large bays is Souza [START_REF] De Souza Machado | Metal fate and effects in estuaries: A review and conceptual model for better understanding of toxicity[END_REF][START_REF] Turner | Suspended Particles: Their Role in Estuarine Biogeochemical Cycles[END_REF], a finer description of Cu speciation/partition and isotope compositions would require additional data characterization of samples (e.g. with synchrotron-based X-ray absorption). Fortunately, previous application of selective chemical leaching (ascorbate, acid dilute and alkaline solution) have already indicated that Fe/Mn oxy-hydroxides and sulphides are the main carrier phases controlling Cu mobility in Toulon bay (Dang et al., 2015). Clay minerals play this important role on Cu dynamics due to their abundance and high charge surface. In these solid phases, Cu may occur within the lattice of aluminosilicates, adsorbed on the surface of clays and amorphous solids or complexed by organic matter [START_REF] Sparks | Toxic Metals in the Environment: The Role of Surfaces[END_REF]. As a trace component of silicates, Cu probably inherits the lithogenic Cu isotope composition close to the average of the Upper Continental Crust (UCC, +0.07± 0.10‰, [START_REF] Liu | Copper isotopic composition of the silicate Earth[END_REF][START_REF] Savage | Copper isotope evidence for large-scale sulphide fractionation during Earth's differentiation[END_REF]. For other carrier phases, laboratory experiments under controlled conditions show a preferential association of the heavy isotope into solid surfaces, such as oxy-hydroxides [START_REF] Balistrieri | Fractionation of Cu and Zn isotopes during adsorption onto amorphous Fe(III) oxyhydroxide: Experimental mixing of acid rock drainage and ambient river water[END_REF][START_REF] Pokrovsky | Copper isotope fractionation during its interaction with soil and aquatic microorganisms and metal oxy(hydr)oxides: Possible structural control[END_REF] and organic complexes [START_REF] Bigalke | Copper Isotope Fractionation during Complexation with Insolubilized Humic Acid[END_REF]. Therefore, the natural Cu partition in sediments and associated isotope signatures of the Large bay seems consistent with previous field and laboratory findings reported in the literature.

On the other hand, the outlier δ 65 Cu NIST values corresponding to the sampling stations 25 and 26 are similar to the light isotope signatures found in the Small bay. This can be related to a zone of dredge spoils near St. Mandrier (Fig. 1) or a preferential hydrodynamic deposition of fine particles enriched in contaminants in this area [START_REF] Tessier | Study of the spatial and historical distribution of sediment inorganic contamination in the Toulon bay (France)[END_REF]. This area was used in the 1980's for contaminated spoils from the north of the Small bay.

Subsequent transport and dispersion promoted by reworking and resuspension (large boat traffic, storm events) of the Cu-enriched dredged sediments could explain the anomalous value of -0.79 at St.26.

The dominant anthropogenic Cu signature in the Small bay (Fig. 4) tending to negative δ 65 Cu NIST values can result from the mixing of ancient and modern Cu anthropogenic sources by reworking of sediments. It potentially includes a broad spectrum of Cu materials employed in harbor and Navy activities, including alloys, plumbing, antifouling paints, warfare, etc.

In turn, the sample at MIS station, collected in the military harbor, and that presents the highest Cu concentration in the sampling grid, constitutes a remarkably punctual anthropogenic source with a positive δ 65 Cu NIST = +0.19 ±0.05‰ (2s). The isotope signature could be attributed to the burning of warfare material enriched in Cu (explosives, driving bands, fuses and shell casings) during WWII. In 1942, The French Navy scuttled more than navy fleet to avoid its capture by the Nazi forces (Grasset, 2011). It is evident that such source attribution requires additional Cu isotope characterization of ancient and modern industrialized materials as well historical information of manufacturing processes and ore purchasing. High contamination of Cu and other metals induced by bullet shells have been reported for battlefields of the Great War: in Verdun (France), burning of arsenic-containing shells resulted in anomalous concentrations in soils for Cu (9113 µg.g -1 ), Pb (5777 µg.g -1 ) and Zn (90190 µg.g -1 , Bausinger et al., 2007), while in the region around Ypres (Belgium) was estimated an input of 2800 t of Cu in the upper 0.5 m of soil due to shelling accumulation (Meirvernne et al., 2008). These findings support our preliminary interpretation of Cu isotope signatures in the military zone of Toulon bay.

Similarly to "MIS" sample, the "Petite Mer" and "St.16" samples present also relatively heavy δ 65 Cu NIST values. However, Petite Mer sample has a very low Cu content (ca.

15 µg.g -1 , EF Cu = 1.9) and a typical radiogenic 206 Pb/ 207 Pb ratio (1.1801), which supports a natural Cu origin, as lattice component of minerals. A possible explanation for its heavy δ 65 Cu NIST and low Cu concentration value would be a granulometry effect associated to coarser mineral grains (e.g, quartz) with low mineral surface charge that sorbs anthropogenic Cu [START_REF] Sparks | Toxic Metals in the Environment: The Role of Surfaces[END_REF]. Instead, the sampling St.16 (c.a, 44, µg.g -1 , EF Cu = 5.0) is moderately enriched, suggesting that its heavy δ 65 Cu NIST value (+0.34‰) reflects an anthropogenic origin.

Several sedimentary archives studied previously by Dang et al. 2015a have showed a remarkable peak of Cu concentration corresponding to WWII events followed by a continuous decreasing of Cu inputs in the period post-WWII. However, the sediment core collected near this zone presented no attenuation of Cu concentrations in the profile. This could be explained by sediment reworking of Cu-enriched layers and mixing of Cu isotope compositions related to past events of WWII impacts. Alternatively, a punctual unknown anthropogenic source cannot be discarded. A temporal analysis of Cu isotopes will be useful to discriminate past and present anthropogenic inputs, and hence, elucidating this question.

Zn isotope systematic

Zinc concentrations varied from 19 to 1471 µg. g -1 , while δ 66 Zn JMC ranged from -0.27 to +0.28‰ (Table 1, Fig. 4). Unlike Cu and Pb isotope systems, Zn isotope compositions in the Small and Large bays greatly overlap: +0.06 ±0.05 (1s) and +0.06 ±0.11(1s), respectively.

Indeed, δ 66 Zn JMC and Zn concentrations do not correlate (Fig. 6a) neither in the Small or Large bay, limiting interpretations in terms of source mixing processes. Interestingly, Zn concentrations correlate very well with Pb concentrations and 206 Pb/ 207 Pb ratios (Fig. 6b,c). [START_REF] Cloquet | Isotopic Composition of Zn and Pb Atmospheric Depositions in an Urban/Periurban Area of Northeastern France[END_REF], or slags (δ 66 Zn > 0.5‰, Yin et al., 2005;Araújo et al., 2017a,b). The final refined Zn metal (Zn 0 ) approximates the average Zn isotope composition of ores [START_REF] John | The isotopic composition of some common forms of anthropogenic zinc[END_REF], when Zn processing yields in metallurgical plants are higher than 95% [START_REF] Sivry | Zn isotopes as tracers of anthropogenic pollution from Zn-ore smelters The Riou Mort-Lot River system[END_REF]Yin et al., 2015). This implies that Zn anthropogenic signatures of alloys used in shipyard or in shells are similar and thus, not discriminable by Zn isotopes [START_REF] John | The isotopic composition of some common forms of anthropogenic zinc[END_REF]Yin et al., 2015).

Natural Zn inputs in coastal systems are associated to inorganic particles issued from weathered soils and lithogenic materials, detrital and authigenic organic matter and in-situ formed minerals, as carbonates and solid sulfide minerals (Araújo et al., 2017a;[START_REF] Mortimer | Metal Speciation (Cu, Zn, Pb, Cd) and Organic Matter in Oxic to Suboxic Salt Marsh Sediments, Severn Estuary, Southwest Britain[END_REF][START_REF] Zhang | Historical records, distributions and sources of mercury and zinc in sediments of East China sea: Implication from stable isotopic compositions[END_REF] [START_REF] Araújo | Assessment of the metal contamination evolution in the Loire estuary using Cu and Zn stable isotopes and geochemical data in sediments[END_REF]. This slight enrichment in light isotopes could be attributed to the regional geology.

Overall, igneous, sedimentary and metamorphic rocks worldwide present homogenous isotope compositions clustering around +0.2 to +0.3‰ [START_REF] Moynier | The Isotope Geochemistry of Zinc and Copper[END_REF], but discrete and detectable Zn isotope heterogeneities on geological settings may occur as result of magmatic differentiation, weathering, crustal subduction processes, mantle melting and volcanisms, and past biological activity [START_REF] Chen | Zinc isotope fractionation during magmatic differentiation and the isotopic composition of the bulk Earth[END_REF][START_REF] Doucet | Zn isotope heterogeneity in the continental lithosphere: New evidence from Archean granitoids of the northern Kaapvaal craton, South Africa[END_REF][START_REF] Inglis | The behavior of iron and zinc stable isotopes accompanying the subduction of mafic oceanic crust: A case study from Western Alpine ophiolites: Fe and Zn isotopes during subduction[END_REF][START_REF] John | Zinc and cadmium stable isotopes in the geological record: A case study from the post-snowball Earth Nuccaleena cap dolostone[END_REF][START_REF] Kunzmann | Zn isotope evidence for immediate resumption of primary productivity after snowball Earth[END_REF][START_REF] Lv | Copper and zinc isotope fractionation during deposition and weathering of highly metalliferous black shales in central China[END_REF][START_REF] Wang | Zinc isotope fractionation during mantle melting and constraints on the Zn isotope composition of Earth's upper mantle[END_REF]. Recent studies have demonstrated how Zn isotope heterogeneity in the continental lithosphere is useful to elucidate petrogenesis evolution, e.g, for discriminating the influence of differentiation processes from melting products in felsic rocks [START_REF] Doucet | Zn isotope heterogeneity in the continental lithosphere: New evidence from Archean granitoids of the northern Kaapvaal craton, South Africa[END_REF]. of these rocks lead to the progressive breakdown of primary minerals and the associated formation of phyllosilicates and iron oxides depleted in Zn content [START_REF] Suhr | Elemental and isotopic behaviour of Zn in Deccan basalt weathering profiles: Chemical weathering from bedrock to laterite and links to Zn deficiency in tropical soils[END_REF].

Changes on Zn-mineral bearing phases combined with biotically-and kineticallycontrolled sorption reactions in Fe-oxyhydroxides and clay minerals surfaces could explain the enrichment of light Zn isotopes present in the weathered residual material [START_REF] Guinoiseau | Zinc and copper behaviour at the soil-river interface: New insights by Zn and Cu isotopes in the organic-rich Rio Negro basin[END_REF][START_REF] Guinoiseau | Zn Isotope Fractionation during Sorption onto Kaolinite[END_REF][START_REF] Little | The oceanic mass balance of copper and zinc isotopes, investigated by analysis of their inputs, and outputs to ferromanganese oxide sediments[END_REF][START_REF] Suhr | Elemental and isotopic behaviour of Zn in Deccan basalt weathering profiles: Chemical weathering from bedrock to laterite and links to Zn deficiency in tropical soils[END_REF]. In turn, the heavy Zn isotope would be leached and enriched in the dissolved phase [START_REF] Little | Key role of continental margin sediments in the oceanic mass balance of Zn and Zn isotopes[END_REF]. Further studies focused in the Mediterranean French geology will be helpful to establish the baseline of Zn isotopes in the region.

Conclusions

Toulon bay surface sediments present a remarkably anthropogenic contamination for Cu, Zn, Pb, Ag and Cd elements. Overall, these anthropogenic elements display a similar spatial distribution in the Toulon sub-bays, with the highest levels of contamination occurring in the Small bay and the lowest in the Large bay. Multivariate and geostatistics approaches show a clear inter-bay differences, highlighting the highest concentrations for Cu, Zn, Pb and Ag in proximity to civil and military ports and the Las river outlet. However, these methods cannot deconvolute the multiple origins of Pb, Zn, and Cu.

The Pb isotope systematics indicate a binary source mixing process involving industrial and natural sources, while Cu isotope records discriminated, at least two distinct anthropogenic signatures, preliminarily interpreted as associated to shipyard alloys and warfare materials. In turn, Zn isotope composition of natural and anthropogenic sources overlap. For Zn, analysis of specifics mineral phases of sediments rather than bulk could be more effective to decipher its natural and anthropogenic origins.

This study presents the first detailed spatial distribution of Pb, Cu and Zn isotope compositions for an aquatic system and demonstrates the use of Cu isotopes as tracers of anthropogenic sources in coastal environments. Considering the crucial role of surface sediments in the transfer of metal contaminants to water column and biota, this first isotope characterization will be useful to track the metal contaminant pathways from abiotic compartments to marine food webs in future studies. Further investigations on Cu, Zn and Pb isotope compositions of surface sediments will be useful to identify potential new anthropogenic inputs, and thus, to better assess the change of source apportionment since WWII.

  associations in the Small and Large bays. Supported by EF values, these associations show that: (i) Cu, Zn, Pb, Ag and Cd represent elements of dominant anthropogenic origin and are grouped along the first component (x-axis); (ii) Li, Co, V, Cr and Ni represent elements predominantly issued from natural sources and are grouped along the second component (yaxis); and (iii) that Mn stands apart in relation to the other elemental groups representing an element of dominant natural origin influenced by diagenetic processes. The diagenesis occurring in sediments under anoxic conditions in the Small bay reduces Mn oxyhydroxides

  noticeable in the Fig. 4. In the plot of δ 65 Cu NIST vs. 1/[Cu] (Fig. 5b) well-defined data clusters are identified for each bay. However, in both sub-bays isotopic outliers are observed: the Small bay clusters overall light Cu isotope compositions, but the sampling stations MIS (δ 65 Cu NIST = +019), Petit Mer (δ 65 Cu NIST = +0.01) and St.16 (δ 65 Cu NIST = +0.34) stand out with positive δ-values (Fig. 5b). In turn, the Large bay is most dominated by positive and narrow range of δ 65 Cu NIST values, from -0.04 to +0.24‰, but the negative δ-values of sampling St.25 (δ 65 Cu NIST = -0.22) and 26 (δ 65 Cu NIST = -0.79) diverge from the Large bay clustering data (Fig.5b). In terms of EF Cu values, the St.37 and 42 differ from the rest of sediments of the Large bay. Unfortunately, we do not dispose of their Cu isotope compositions. Disregarding these outlier data, a general trend toward lighter δ 65 Cu NIST associated with an increase of Cu concentrations or EF Cu (Fig.5c) is noticeable, which suggests a source mixing process involving anthropogenic and natural sources.The overall positive and heavier δ 65 Cu NIST values and associated low Cu EF values in the Large bay reflect a dominant signature of natural Cu. As coastal sediments integrate inorganic and organic compounds derived from detrital and autochthonous materialssuch as, clay minerals, amorphous solids, biogenic carbonates, organic debris-(Bianchi, 2007; de

  The local geology of Toulon bay basin encompasses mainly Permian and Triassic sediments, Paleozoic metamorphic rocks (www.geoportail.gouv.fr). Evoking possible mechanism of Zn isotope fractionation during high-temperature such as magmatic differentiation could explain this slight tendency to lighter isotope values in the local geology of Toulon bay. Alternatively, this light Zn isotope geochemical baseline of sediments of Toulon bay could be associated to the contribution of sedimentary rocks sources. Weathering processes occurring in the formation

Fig. 1 .

 1 Fig. 1. Map of Toulon bay and sampling sites of surface sediments.

Fig. 2 .

 2 Fig. 2. Interpolation maps of trace elements enrichment factors (EF) in surface sediments. Enrichment factor values were classed as suggested in other studies (Chen et al., 2007): EF < 1 indicates no enrichment; EF < 3 for minor enrichment; EF = 3-5 for moderate enrichment; EF = 5-10 for moderately severe enrichment; EF = 10-25 for severe enrichment; EF = 25-50 for very severe enrichment; and EF > 50 for extremely severe enrichment.

Fig. 3 .

 3 Fig. 3. Principal component analysis of surface sediments from Toulon bay: (a) the main components and factor loads; (b) the sample factor scores.

Fig. 4 .

 4 Fig. 4. Interpolation map of isotope compositions of Pb, Cu and Zn in surface sediments.

Fig. 5 .

 5 Fig. 5. Plots of elemental and isotope data of Cu and Pb. Copper concentrations [Cu] are expressed in µg. g -1 .

Fig. 6 .

 6 Fig. 6. Plots of elemental and isotope data of Zn. Zinc concentrations [Zn] are expressed in µg. g -1 .

  This allows using Pb isotopes to identify the best representative Zn anthropogenic and Zn natural samples. Using the Pb isotope ratios associated to natural ( 206 Pb/ 207 Pb ~ 1.180) and anthropogenic ( 206 Pb/ 207 Pb ~ 1.165) sources as cut-off lines (Fig.6d), we sort them in two groups. The natural Zn samples with Zn concentrations lower than 60 µg.g -1 and a δ 66 Zn JMC

	average of +0.14 ±012 ‰ (2s, n= 9), and anthropogenic Zn samples with average
	concentrations of 383 µg.g -1 and δ 66 Zn JMC values of +0.04±0.08 ‰ (2s, n= 12). These
	δ 66 Zn JMC values estimated as the best representative of natural and anthropogenic Zn sources
	are statically equal (test t, p >0.5) and show an overlap between isotope signatures of natural
	and anthropogenic sources in Toulon bay.
	It is noted that the estimated anthropogenic δ 66 Zn JMC values of +0.04 ±0.08 (2s, n= 12)
	are close to the isotope range reported for Zn sulfide (from +0.13 to +0.15‰) and silicates
	ores (-0.06 to +0.41‰, (Mondillo et al., 2018), and are clearly distinguished from those
	attributed to industrial atmospheric emissions (< -0.6‰, Souto-Oliveira et al., 2018; from -
	0.52 to +0.02‰, Mattielli et al., 2009; from -0.27 to +0.15‰, Ochoa et al., 2016; from
	+0.16 to +0.63‰,

  . The detrital-inorganic materials dominate sediment compositions of Toulon bay. However, the estimated natural δ 66 Zn JMC values of +0.14 ±0.12‰ are slightly lighter than the Zn isotopes values reported to the Upper Earth Crust,

	estimated as +0.28 ±0.05‰ (Chen et al., 2013) and other unpolluted sediments from estuaries
	worldwide: (+0.28‰, Sepetiba bay, Araújo et al., 2017a; +0.3‰, Gironde estuary, Petit et al.,
	2015 ; +0.3‰, the Loire estuary,

Table 1 .

 1 Copper, zinc and lead data description of Toulon bay sediments.

				All Toulon bay
		n	Range	Minimum Maximum	Mean	Std. Deviation (1s)
	δ 65 Cu (‰)	44	1.13	-0.79	0.34	-0.05	0.24
	δ 66 Zn (‰)	52	0.55	-0.27	0.28	0.06	0.09
	206 Pb/ 207 Pb	52	0.0407	1.1591	1.1998	1.1717	0.0081
	Cu (µg.g -1 )	52	662	4	666	85	111
	Zn (µg.g -1 )	52	1452	19	1471	178	223
	Pb (µg.g -1 )	52	456	13	470	93	83
				Small bay
		n	Range	Minimum Maximum	Mean	Std. Deviation (1s)
	δ 65 Cu (‰)	23	0.83	-0.49	0.34	-0.17	0.19
	δ 66 Zn (‰)	24	0.18	-0.02	0.16	0.06	0.05
	206 Pb/ 207 Pb	24	0.0210	1.1591	1.1801	1.1664	0.0043
	Cu (µg.g -1 )	24	651	15	666	147	137
	Zn (µg.g -1 )	24	1411	60	1471	302	278
	Pb (µg.g -1 )	24	438	32	470	144	94
				Large bay
		n	Range	Minimum Maximum	Mean	Std. Deviation (1s)
	δ 65 Cu (‰)	21	1.03	-0.79	0.24	0.08	0.23
	δ 66 Zn (‰)	28	0.55	-0.27	0.28	0.06	0.11
	206 Pb/ 207 Pb	28	0.0357	1.1641	1.1998	1.1763	0.0079
	Cu (µg.g -1 )	28	141	4	145	32	32
	Zn (µg.g -1 )	28	214	19	233	72	50
	Pb (µg.g -1 )	28	131	13	145	49	35

Acknowledgements

This work was conducted within the framework of the POLLUSOLS and SCOTTTI projects.

They were funded by the regional government of the Pays de la Loire (France) and AERMC (Agence de l'Eau Rhône Méditerranée Corse), respectively. A part of this research is also a contribution to the MERMEX-Merite research program. This study benefits from sediment sample collection occurring from the CARTOCHIM project funded by "Toulon-Provence-Méditerranée (TPM)", and Agence de L'Eau Rhône Méditerranée Corse. The authors thank the anonymous reviewers for their helpful suggestions and comments on earlier drafts of the manuscript. references Ahlf, W., Drost, W., Heise, S., 2009. Incorporation of metal bioavailability into regulatory frameworksmetal exposure in water and sediment. J Soils Sediments 9, 411-419. https://doi.org/10.1007/s11368-009-0109-6 Araújo, D., Machado, W., Weiss, D., Mulholland, D.S., Boaventura, G.R., Viers, J., Garnier, J., Dantas, E.L., Babinski, M., 2017b. A critical examination of the possible application of zinc stable isotope ratios in bivalve mollusks and suspended particulate matter to trace zinc pollution in a tropical estuary. Environmental Pollution 226, 41-47. https://doi.org/10.1016/j.envpol.2017.04.011 Araújo, Daniel F., Boaventura, G.R., Machado, W., Viers, J., Weiss, D., Patchineelam, S.R., Ruiz, I., Rodrigues, A.P.C., Babinski, M., Dantas, E., 2017a. Tracing of anthropogenic zinc sources in coastal environments using stable isotope composition. Chemical Geology 449, 226-235. https://doi.org/10.1016/j.chemgeo.2016.12.004 Araújo, D.F., Boaventura, G.R., Viers, J., Mulholland, D.S., Weiss, D., Araújo, D., Lima, B., Ruiz, I., Machado, W., Babinski, M., Dantas, E., 2016. Ion Exchange Chromatography and Mass Bias Correction for Accurate and Precise Zn Isotope Ratio Measurements in Environmental Reference Materials by MC-ICP-MS. Journal of the Brazilian Chemical Society. https://doi.org/10.5935/0103-5053.20160167