Charles Peureux

Fabrice Ardhuin
email: fabriceardhuin@ifremer.fr

Pedro Veras Guimarães

On the Unsteady Steepening of Short Gravity Waves Near the Crests of Longer Waves in the Absence of Generation or Dissipation

Keywords:

crest of the longer waves

Introduction

The evolution of random wave fields is a complex problem that has often been treated as a spatially homogeneous or stationary process in which non-linear wave evolution is represented by 4-wave nonlinear interactions that control the cascade of energy to shorter components with a source of energy from the wind and dissipation by wave breaking (e.g., [START_REF] Hasselmann | On the non-linear energy transfer in a gravity wave spectrum, part 1: general theory[END_REF][START_REF] Phillips | Spectral and statistical properties of the equilibrium range in wind-generated gravity waves[END_REF]. A particular difficulty of the problem is the interaction of multiple scales that introduce effective inhomogeneities, in particular the short waves are expected to respond to the the orbital velocities and apparent gravity of the long waves. Recent parameterizations of the spectral wave dissipation have used empirical representations of enhanced breaking of short waves induced by long waves [START_REF] Donelan | Modeling waves and wind stress[END_REF][START_REF] Romero | Distribution of surface wave breaking fronts[END_REF]. This enhanced breaking may explain a large change in the energy level of waves with periods 1 to 2 s. However, it requires a detailed justification. Here we will ignore the effect of the short wave non-linearity that are discussed by [START_REF] Zhang | On the stability of weakly nonlinear short waves on finite-amplitude long gravity waves[END_REF], and that are certainly important for the eventual evolution of short waves toward breaking. Instead we focus on the non-homogeneity of linear shorter wave components that is caused by the presence of linear longer wind-generated waves. This effect is generally called modulation. Very similar effects occur when surface gravity waves interact with internal waves [START_REF] Lewis | On the interaction of internal waves and surface gravity waves[END_REF] or much longer gravity waves such as tides [START_REF] Ardhuin | Numerical wave modeling in conditions with strong currents: dissipation, refraction and relative wind[END_REF].

Much of the the previous modulation theory and observation work was motivated by the analysis of remote sensing data, in particular the modulated power of microwave radar signals backscattered from short gravity-capillary waves, which is routinely used to measure wind and waves. These investigations have produced very useful theoretical and empirical transfer functions in which the local wind speed and direction play an important role (e.g., [START_REF] Keller | Microwave scattering and the straining of wind-generated waves[END_REF][START_REF] Hasselmann | Theory of Synthetic Aperture Radar ocean imaging: a MARSEN view[END_REF][START_REF] Hara | Observation of hydrodynamic modulation of gravity-capillary waves by dominant gravity waves[END_REF]. The early works were also associated to possible net growth or decay of long waves as an indirect interaction of long waves and wind through the presence of short waves [START_REF] Longuet-Higgins | Action of a variable stress at the surface of water waves[END_REF][START_REF] Hasselmann | On the mass and momentum transfer between short gravity waves and larger-scale motions[END_REF][START_REF] Garrett | On the interaction between long and short surface waves[END_REF]. Such effects may indeed be relevant to swell dissipation or the generation of the long-period fore-runners described by [START_REF] Munk | Directional recording of swell from distant storms[END_REF].

Our purpose differs from these previous works as we focus on not-so-short modulated components, with typical wavelengths around 1 m. For these wavelengths, and for the average wind speed of 7 m/s, a typical growth rate given by eq. ( 1) in [START_REF] Plant | A relationship between wind stress and wave slope[END_REF] is β ≈ 0.01 s -1 . Over a typical long wave period of 10 s, this gives a growth of only 10%, which we shall neglect. Our interest is the intermittent variability of wave properties, including the mean energy level of the not-so-short waves as a function of direction. One particular motivation for the present work is that the energy balance of waves beyond 2 times the dominant frequency requires some poorly understood "cumulative dissipation" (e.g., [START_REF] Banner | Multiscale measurement of ocean wave breaking probability[END_REF][START_REF] Banner | Refined source terms in wind wave models with explicit wave breaking prediction. part I: Model framework and validation against field data[END_REF][START_REF] Ardhuin | Semi-empirical dissipation source functions for wind-wave models: part I, definition, calibration and validation[END_REF], which may be related to modulation-type processes, possibly associated with steep long waves [START_REF] Kharif | Some aspects of the kinematics of short waves over longer gravity waves on deep water[END_REF]. Recent observations of meter-scale waves have revealed a suppressed energy level for short waves aligned with the dominant waves [START_REF] Leckler | Analysis and interpretation of frequency-wavenumber spectra of young wind waves[END_REF][START_REF] Peureux | Note on the directional properties of meter-scale gravity waves[END_REF], that cannot be explained by nonlinear 4-wave interactions alone [START_REF] Peureux | Observation et modélisation des propriétés directionnelles des ondes de gravité courtes[END_REF].

It has been recognized for a long time that the dominant waves may organize the evolution of the shorter components by many modulation processes (Longuet-Higgins & Stewart, 1960, hereinafter LH&S), with a possible influence on the spectra of both short and long waves [START_REF] Hasselmann | On the mass and momentum transfer between short gravity waves and larger-scale motions[END_REF][START_REF] Banner | Wavenumber spectra of short gravity waves[END_REF]. We have thus decided to revisit the modulation effect for meter-scale short gravity waves in the presence of longer waves. Here we consider a highly simplified model without sources or sinks of energy, following the classic paper of LH&S. This conservative model is known to be unrealistic in the case of waves over tidal currents, for which wave breaking severely limits the predicted steepness [START_REF] Ardhuin | Numerical wave modeling in conditions with strong currents: dissipation, refraction and relative wind[END_REF]. Still it exhibits some interesting behavior that has not been described before. In particular we find that the short wave steepness can be much larger than predicted by LH&S.

The theoretical model is described in section 2. Section 3 presents the solutions and describes the importance of the terms neglected by LH&S. Section 4 gives a prediction for an envelope instability that is verified numerically. These results are then discussed in section 5 in the context of random wave breaking and a cumulative effect of long waves on short waves.

Problem Statement

Here we consider a narrow-band field of short waves propagating over a varying large scale velocity field with a horizontal component U and a vertical component W in deep water. The short waves are characterized by an energy density E per unit horizontal surface, a wave vector k, an intrinsic radian frequency σ = σ(k), and an associated action density (e.g. [START_REF] Phillips | The dynamics of the upper ocean[END_REF])

N = E σ . (1)
The short waves have intrinsic phase velocity C = σ/k and group velocity C g = ∂σ/∂k. The distribution of the short wave energy in space and time is given by the waveaction conservation equation and the conservation of crests [START_REF] Bretherton | Wavetrains in inhomogeneous moving media[END_REF][START_REF] Phillips | The dynamics of the upper ocean[END_REF],

∂ t N + ∇ • [(C g + U ) N ] = 0 (2) ∂ t k + ∇ (σ + k • U ) = 0, (3)
where ∇ = [∂ x , ∂ y ], provided the external current varies slowly compared to the short waves wave length, and neglecting other sources and sinks of energy such as wind input or breaking dissipation. The short wave energy is related to the local envelope amplitude a as given by eq. (4.11) in LH&S,

E = 1 2 ρga 2 1 + 1 2g ∂ t W . (4)
Here we assume that the short waves are in deep water, so that the radian frequency σ and wave number k are related by the dispersion relation

σ 2 = g k, (5)
where the apparent gravity g is the gravity acceleration modified by the long wave induced downward vertical acceleration

g = g 1 + 1 g ∂ t W . (6)
which slightly differs from the factor in equation ( 4). This was not taken into account in LH&S but it is taken into account in later works (e.g. [START_REF] Longuet-Higgins | A stochastic model of sea-surface roughness[END_REF][START_REF] Elfouhaily | Improved electromagnetic bias theory: Inclusion of hydrodynamic modulations[END_REF]. In deep water, the group speed is related to the phase velocity through

C g = C 2 = 1 2 g k . ( 7
)
Together with ( 5), ( 6) and (7), the system of equations ( 2)-(3) governs the evolution of the short waves envelope.

We consider a large-scale flow U given by a deep-water monochromatic longer wave of steepness ε L 1, wave number k L , radian frequency σ L , phase speed

C L = σ L /k L and a phase ψ = k L x -σ L t U = ε L C L cos ψ, (8)
the vertical velocity

W = ε L C L sin ψ (9)
and the apparent gravity

g = g (1 + ε L M g cos ψ) , (10)
where M g = -1 [START_REF] Phillips | The dynamics of the upper ocean[END_REF]. The long wave orbital velocity given by eq. ( 8) is associated with a long wave surface elevation ε L /k L cos ψ, so that ψ = 0 • corresponds to the crests and ψ = 180 • to the troughs.

LH&S proposed solutions for the 1D case in which they neglected the gradients of the short waves group velocity compared to the gradients of the orbital currents. It simplified greatly the algebra and lead them to a general first order solution for any short waves modulated variable

X(ψ) = X i (1 + ε L M X cos ψ) , (11)
where the property X of the short waves only depends on the long wave phase and X i is the mean value of X. M X is usually called the modulation transfer function (hereinafter MTF) for the variable X. M X is of the order of unity for most variables, especially M N = 1 and M k = 1 in the LH&S approximate solution.

In the following, we investigate the modulation in more general conditions, and this departure from the mean short waves behavior, normalized by the long wave steepness, is a periodic function of time and space,

∆X(x, t) = X(x, t) -X i ε L X i , (12)
and the MTF M X in (11) appears in the Fourier series expansion of ∆X as a ψ-periodic function:

∆X(x, t) = M X cos (ψ + ϕ X ) + ∞ n=2 M X,n cos [nψ + ϕ X,n ] (13)
This generalization of the modulation is particularly necessary when the modulation is not stationary. From now on, we will often refer to the maximum modulation after a given time to characterize the short waves field.

For simplicity, we consider the problem in which all quantities are invariant along the y-axis, i.e. along the long wave crests, so that k y remains constant along the propagation, and all y partial derivatives are cancelled. This assumption transforms equations ( 2) and (3) to their simplified versions

∂ t N + ∂ x [(γC/2 + U ) N ] = 0 (14) ∂ t k x + ∂ x [(C/γ + U ) k x ] = 0, (15)
where γ = k x /k is the cosine of the short waves propagation direction, which will be used for numerical simulations. Changing variables from N and k x to N and C, with C = g /k, gives the alternative set of equations, which will be used for analytical calculations,

∂ t + γ 2 C + U ∂ x N = -N (3γ 2 -2)/(2γ)∂ x C -N (1 -γ 2 )C/(2γg )∂ x g -N ∂ x U (16) ∂ t + γ 2 C + U ∂ x C = (γ 2 C/2)∂ x U + C/(2g ) [∂ t + (γC + U ) ∂ x ] g , (17)
with derivation details given in appendix A.

At time t = 0, a monochromatic and unidirectional short wave is considered, with a slight deviation along the long wave direction that can encompass any kind of random modulation by the latter, considering the modulation of the order of the long wave steepness, in other words

N (t = 0, x) = N i + ε L δN (x) (18)
C(t = 0, x) = C i + ε L δC(x) (19)
k x (t = 0, x) = γ i k i + ε L δk x (x) (20) with δC(x) = 1 2 M g cos ψ -γ i δk x (x) k i + O ε 2 L . ( 21
)
We may also consider more general initial conditions,

δN (t = 0) /N i = M N,i cos ψ + ∞ n=2 (M N ) n,i cos nψ (22) δk x (t = 0) = γ i k i M k,i cos ψ + ∞ n=2 (M k ) n,i cos nψ . ( 23
) LH&S initial condition corresponds to M N,i = M k,i = 1 and (M N )n, i = (M k )n, i = 0, n ≥ 2.
In the following, initial conditions harmonics with n ≥ 2 will be set to zero for simplicity: we have found that these give perturbations that do not qualitively change the nature of the solutions.

As a result the initial condition on C is (see appendix A for more details),

C (t = 0) = C i (1 + ε L M C,i cos ψ) (24) at order ε L , where M C,i = 0.5(M g -γ 2 i M k,i ).
To summarize, our problem is fully determined by the initial conditions for the short wave wavenumber and action M k,i (t = 0, x) and M N,i (t = 0, x). and 3 parameters, namely the long wave steepness ε L , the initial scale ratio k i /k L , and the short waves propagation direction with respect to the long wave θ i , or γ i = cos θ i . Simulation results are linear in N i , and independent of the short waves initial steepness.

In the following, the coupled equations ( 14) and ( 15), together with initial conditions ( 22) and ( 23) are solved numerically using a procedure detailed in appendix B and summarized here. As spatially periodic solutions are sought, the simulation spatial domain is reduced to a single long wave wave length, from ψ = 0 • to ψ = 360 • , with periodic boundary conditions. The problem is solved in a frame moving with the long wave phase velocity. Following up on our interest in short wave breaking, we will particularly focus on the long wave crests, ψ = 0 • , because the largest short wave steepness is obtained there.

Modulation of Two Collinear Wave Trains

Numerical Simulation

We first consider the simplified 1D configuration where short and long waves are aligned along the x-axis. The short waves do not refract and γ = γ i = 1 is a constant for the moment. Starting from a uniform distribution, the evolution of short wave action N , wavenumber k x and short waves steepness ak is shown on figure 1 for

k i /k L = 10 and ε L = 0.1.
The space and time evolution of the wavenumber k x is nearly periodic, with higher values (shorter wavelengths) around the crests of the long waves. Away from the long wave crests, there is a local maximum that travels at the speed of a short wave group with a position shown with the dotted line. The pattern of the action and steepness are similar except that they are clearly not periodic: the wave packet that starts from the trough gains action and steepness each time it goes over a long wave crest, but it only loses part of that gain as it goes in the following trough. Conversely, a wave packet that starts from the wave crest loses action and steepness as it goes over a trough, but only recovers some of that when going over the next crest. What happens is a progressive condensation of action along the trajectory of the first short wave group, the one starting from the trough, and the passage of the wave crests. As the total action is conserved, this gain of one wave group goes with a loss of the neighboring wave groups.

Because the long wave orbital contribution converges on the forward face of the long waves, the increase in short wave action occurs one quarter of period later, i.e. on the crest. This effect is already present in the LH&S solution. What we find here, is that the additional change in advection velocity due to the modulation of the short wave group speed gives an additional amplification. This amplification has a non-zero mean effect because the modulations of N and C g = 0.5 g/k are correlated: the average perturbation in the advection of action is a product C g N that on average is not zero because of the phase relation between N and k. As a result every long wave cycle brings a little more convergence and the maximum value of N keeps growing at each cycle. This leads to a pattern of very sharp action maxima that travels along the long wave profile at a speed that is very close to the short wave group speed (dotted lines in figure 1). Every time the maximum goes over the crest of a long wave it is amplified a little more A more quantitative understanding is provided by the time evolution of both the wave action and the wave number modulations as presented on figure 2 following a long wave crest, i.e. ψ = 0 • , corresponding to a vertical slice through figure 1. Two simulations are shown starting from either the modulated LH&S profile (M N,i = M k,i = 1) or a uniform distribution along the long wave profile (M k,i = M N,i = 0). This second option corresponds to the sudden appearance of a long wave perturbation in the middle of a homogeneous short wave field. It is less clear under which the first option may be realized, given that the solution appears to be unstable. Indeed, the top panel suggests that the LH&S solution is much more stable with modulations locked in phase with the long wave crests and an amplitude of 1, but still the numerical solution slowly deviates and produces higher and lower values of N and k.

In general the modulation of k is less pronounced than the modulation of N . The time evolution of k is given by eq. ( 15) alone, as C is a function of k only, and develops into the typical saw-tooth profile due to the non-linearity brought by the piece of the advection velocity (the short wave phase speed) that is proportional to k -0.5 . This is similar to the formation of shocks for shoaling waves in the shallow water equation. Such a formation of shocks is a common property of all nonlinear Burgers-type advection equations [START_REF] Whitham | Linear and nonlinear waves[END_REF]. The numerical scheme here does not fully capture the shocks and keeps them relatively smooth.

The combined behavior of N and k leads to a progressive steepening of the short waves, due to the growth of the wave action peaks almost in phase with those of the wave number (see figure 2). Again, while the steepness of the LH&S initial conditions

[X(¢ = 0 ∘ , t ) X i ] /(ε L X i )
Generalized modulation almost does not evolve, the maximum steepness from any other cosine-distributed initial condition is much larger.

The longer term evolution shown in figure 2 later has M N exceeding 10 after 10 long wave periods. This focusing is typical of Burgers-type advection equations. The later stabilization at M N ≈ 40 for larger times (see fig. 3) is beyond the scope of the present paper and may be due to numerical diffusion. Anyway, such large values should never be reached for short waves with average steepness ak > 0.1 in the presence of long waves with similar steepness, because short wave breaking should occur before.

The growth of the wave action disappears when the advection is linearized, as done by LH&S, when neglecting the ∂ x C g terms in both equations ( 2

Stability Analysis

The behavior of the short waves is investigated analytically in a weak nonlinear interactions framework by switching to the alternative set of equations ( 16) and ( 17), with initial conditions ( 22) and ( 24). The procedure used to solve this problem is detailed in appendix C for the 2D and summarized here in the 1D case. Using a classical perturbation expansion in powers of the long wave steepness ε L , the perturbation induced by the long waves modulation through the first order solution {N 1 , C 1 } on the initial state {N 0 , C 0 } follows (see appendix C)

[∂ t + (C 0 /2)∂ x ] N 1 = -N 0 ∂ x U -(N 0 /2)∂ x C 1 (25) [∂ t + (C 0 /2)∂ x ] C 1 = (C 0 /2)∂ x U + C 0 /(2g) (∂ t + C 0 ∂ x ) g . ( 26
)
Already the sources of the instability are visible. There is a single way coupling occurring here. The evolution of C 1 given by ( 26) is stand-alone, and is plugged into the one of N 1 through the right-hand side of equation ( 25). As both equations have the same left-hand side operator, they have the same natural frequency. The forcing operated by the term ∂ x C 1 in the equation ( 25) forces N 1 at its resonance frequency (secular term), giving rise to an exponential growth of the solution. That term was precisely neglected by LH&S.

In that case, the perturbative expansion cannot be applied, and approximate solutions are obtained using a multiple time scale method (see appendix C). Especially, the resonance results in an exponential growth of the wave action maxima (the total wave action being conserved) with typical e-folding time τ estimated from (C28), which reduces to

1 τ = |Ω| 1 4 ε L σ L C i C L |1 -M k,i | (27)
assuming a 1D case and using C i /C L 1 . This expression is valid for times O(ε -1 L ), typically of the order 10 long wave periods. The growth which takes place does not

-2 -1 0 1 2 3 4 M N, i -2 -1 0 1 2 3 4 M k, i Numerical model -2 -1 0 1 2 3 4 M N, i
S abili y analysis change the total wave action over the long wave profile, but modifies its distribution as explained in the previous section. A comparison of the approximate solution with the numerical simulation of the previous paragraph is presented on figure 3 (gray lines).

For most cases, the growth predicted theoretically is first slower than the one modelled, overpasses it at large enough times, but still presents a behavior to the one computed numerically.

The influence of initial conditions on the growth of the solution is presented in figure 4. For both the numerical and the stability analysis solution, this diagram clearly highlights an unstable equilibrium position around the LH&S initial conditions. More precisely, for the analytical solution, a stationary solution is obtained when

M k,i = M g -2 µ χ 1, (28)
and the oscillations of the action profile are cancelled when

M N,i = µ χ 1, (29)
where µ, χ and µ are defined in appendix C, equations (C25), (C26) and (C30). The combination of ( 28) and ( 29) is very close to the LH&S solution. The exact stationary profile in 1D was studied by [START_REF] Kharif | Some aspects of the kinematics of short waves over longer gravity waves on deep water[END_REF] and deviates slightly from a linearly modulated profile such as the ones taken as initial conditions. From figure 4, the furthest the initial conditions from the stationary solution, the larger the oscillations of the solution around the stationary profile, and the largest the growth of action. It has to be noted that the modulation on k has more significant consequences on the growth of action than the modulation of action itself.

Modulation of Non-Collinear Wave Trains

Numerical simulations and analytical calculations, detailed in Appendix C, have been performed for short waves propagating at an angle relative to the long wave 0 50 100 150 propagation, i.e. with γ < 1. Both results are illustrated in figure 5. The calculations predict a growth of the wave action at the rate

θ i [degrees] 2 4 6 8 10 max{[N(ψ, t < 10T L ) -N i ]/(ε L N i )} ε L = 5% ε L = 10%
1 τ = |Ω| 1 4 ε L σ L C i C L γ i (1 -M k,i ) 3γ 2 i -2 , ( 30
)
as given by eq. (C28), with

C i /C L 1.
The maxima of wave action modulation after 10 long wave periods are presented for two different long wave steepnesses ε L = 0.05 and ε L = 0.1 (solid lines), and compared to the predictions of our analytical model, equation (C1) (dashed lines). The same qualitative behaviors are observed for both the analytical and the numerical solutions. The largest errors of the analytical model are found where the growth rate is the largest. Indeed, the exponential growth expected in the stability analysis generally underestimates the growth of the numerical solution at early times. However, the general directional distribution is well captured. Especially, the growth rate (30) becomes zero at γ i = ± 2/3 which corresponds to θ i 35 • and θ i 145 • . These angles define a ±35 • sector around the upwave and downwave directions for which an excess breaking probability of short waves can be expected.

Qualitatively, the behavior of the non-collinear solutions is similar to collinear wave trains. There still exists for each initial direction θ i a stationary solution, for which modulation does not change, given at order ε L by expressions ( 28) and ( 29). For given cosine-distributed initial conditions, the distance of the initial modulation from the stationary solution determines both the amplitude of the subsequent modulation and instability.

Discussion and Conclusions

From our analysis we found that short waves modulation by monochromatic long waves can be more complex than the steady cosine shape variation in wave action and wavenumber given by the linearized solution of [START_REF] Longuet-Higgins | Changes in the form of short gravity waves on long waves and tidal currents[END_REF]. For a wide class of initial conditions, the short wave energy condenses into very steep and narrow wave packets around the long wave crests. This resonance mechanism is expected to occur for other kinds of dispersive wave interactions, in particular for gravity waves in intermediate water depth. It can also be expected for any kind of dispersive wave packet periodically forced in space and time by an external source. A similar mechanism has already been reported and measured experimentally for the interaction between a surface and an internal gravity waves by [START_REF] Lewis | On the interaction of internal waves and surface gravity waves[END_REF], with a maximum of the interaction when the phase speed of the long wave matches the group speed of the short waves, corresponding to χ = 0 in the analysis reported in Appendix C.

Although the qualitative behavior is the same regardless of the initial conditions, both the oscillations and the growth rate (27) depend on the initial conditions, especially on the magnitude of the short wave perturbation along the long wave profile. There exists a stationary profile of short waves, close to the LH&S solution, but any deviation from that profile results in an exponential growth of the wave action maxima. Ultimately, this should lead to a steepening of the short waves profile and an increase of the breaking probability.

For a wide range of initial conditions there is locally, after about 10 long waves periods, an effect of the long waves on the short waves that is 10 times larger than expected by LH&S. This effect is maximum for short and long wave propagating in the same direction, and is greatly reduced for direction differences larger than 20 degrees. Over this time scale, considering the empirical growth term of [START_REF] Plant | A relationship between wind stress and wave slope[END_REF], the windwave generation can only be neglected for winds under 3 m/s or short wavelengths longer than 1 m. In these conditions, the occurrence of extreme short wave amplitude may be governed by the process described here. In order to be applied to the analysis of meter-scale wave breaking in the presence of longer waves, the instability of short wave packets that we describe here should be investigated in a more realistic setting. In particular we expect that replacing the monochromatic long waves by random long wave groups will limit the magnitude of the maximum steepness as the forcing will only be coherent over only 10 or so long wavelengths. Also, this model may not be valid for all waves in a broad wind wave spectrum, as a scale separation assumption is required. Finally, the present theoretical framework is limited to order ε 1 L solutions, and cannot be applied at long time scales for which nonlinear 4-wave resonant interactions will likely dominate the evolution of the wave field [START_REF] Hasselmann | On the non-linear energy transfer in a gravity wave spectrum, part 1: general theory[END_REF].

In the ocean, short gravity waves ride on a random superposition of waves that is different from the monochromatic waves investigated here. We expect that signif-icant modulations can be generated over the time scale of a few long wave periods, with oscillations around the LH&S solution. Future analysis may further combine generation and dissipation effects with the conservative case described here, possibly following [START_REF] Longuet-Higgins | A stochastic model of sea-surface roughness[END_REF]. The short wave response to wave groups is another topic that should be investigated as a step towards a practical application for the parameterizations of spectral wave evolution. are the maximum wave speeds at a given time step (eigenvalues of the jacobian matrix g ) and where κ R j is computed according to (B9).

The numerical resolution uses 128 points over the long wave wavelength, with periodic boundary conditions q(ψ = 360 • , T ) = q(ψ = 0 • , T ). The time step is chosen so that the Courant number Co is lower or equal to one ∆T = Co∆ψ λ max , (B17) thus insuring the stability of the numerical scheme. As λ max 1, the typical time step is ∆T 1/128 long wave periods.

B2 Test Case

The numerical scheme was tested against the closest analytically solvable problem

∂ t X + ∂ x [U (x, t)X] = 0. (B18)
verified both by X = N and X = k x , which are the same as ( 14) and ( 15) without the self-advection term and the apparent gravity modification, or equivalent to reducing f in (B4) to f (q, ψ) = (-1 + ε L cos ψ) q.

(B19)

With initial condition X(x, t = 0) = X i (1 + ε L M X,i cos ψ) and in the case of a monochromatic long wave orbital current (8), the solution of (B18) is

X (x, t) = X i v [σ L t + u (ψ)] 1 -ε L cos ψ , (B20)
where

u (ψ) = 2/ 1 -ε 2 L arctan (1 + ε L ) / (1 -ε L ) tan (ψ/2) (B21) v (ψ) = 1 + ε L (M X,i -1) r (ψ) -ε 2 L M X,i r 2 (ψ) (B22) r(ψ) = 1 -w 2 (ψ) / 1 + w 2 (ψ) (B23) w (ψ) = (1 -ε L ) / (1 + ε L ) tan 1 -ε 2 L ψ/2 . (B24) (B25)
The classical Lax-Wendroff scheme for nonlinear advection equations [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] has been tested but unrealistically overestimates the growth of instability, contrary to the present scheme. As a drawback, when compared to the analytical solution (B20), the present scheme exhibits a slightly lower accuracy than the Lax-Wendroff scheme. Although conservative, the distributions of N and k x flatten spatially due to MUSCL-RK4. The mean error typically increases by 3% every 100 long waves periods for Co = 1 and 128 spatial samples.

where

µ = 1 - 1 4 C i χC L γ i 3γ 2 i -2 -M g γ i - C i C L 2γ 2 i -1 (C30)
The 1D solutions of section 3 are obtained by setting γ i = 1.

Figure 1 .

 1 Figure 1. Evolution of short wave properties, action N , wavenumber k and slope ak, in the frame of reference moving with the long wave phase ψ, with ψ = 0 corresponding to the long wave crest. The short wave quantities are normalized by the long wave slope εL, following the definition of the modulation in eq. (12). The case shown here is for ki/kL = 10, εL = 0.1, γi = 1 and a Courant number Co=0.1. starting with a uniform distribution of N and k at t = 0. The black dots correspond to successive positions of a linear wave group starting from the trough of the long wave profile and propagating at apparent speed of the short wave groups 0.5Ci (see Appendix C).

Figure 2 .

 2 Figure 2. Time evolution of the relative change of the modulation of wave action N , wave number k and steepness ak following long wave crests for a short waves train with ki/kL = 10, εL = 0.1, γi = 1 and Co = 0.1. Panel (a) corresponds to LH&S initial conditions MN,i = M k,i = 1, and panel (b) to homogeneous initial conditions MN,i = M k,i = 0 as in figure 1. The stability analysis solution and its envelope correspond to the analytical solution of appendix C.

  ) and (3) compared to ∂ x U .

Figure 3 .

 3 Figure 3. Wave action modulation on wave crests (same case as figure 2b for a longer time range).

Figure 4 .

 4 Figure 4. Maximum action modulation after 10 long wave periods from initiation as a function of various initial modulations. The runs are performed with ki/kL = 10, εL = 0.1, γi = 1 and Co = 1. The white crosses stand for LH&S initial conditions MN,i = M k,i = 1.

Figure 5 .

 5 Figure 5. Variation action modulation maxima at the long waves crests over 10 long wave periods, as a function of the difference in direction between the short and long waves. Initial conditions use MN,i = M k,i = 0, εL = 0.1, ki/kL = 10, and we used Courant numbers Co < 0.1. The dashed lines are given by our analytical model.

Acknowledgments

This paper is a purely theoretical work and does not use any data. C.P., F.A. and P.G. were supported by CNES as part of the SKIM preparation program and ANR grants for ISblue (ANR-17-EURE-0015) LabexMER (ANR-10-LABX-19), and MIMOSA (ANR-14-CE01-0012).

Appendix A Alternative Set of Equations A1 Equation for the phase speed C Equation ( 15) is rewritten as a function of C = g /k. From the chain rule, the derivatives of C and k x are connected through dC = ∂C ∂k x dk x + ∂C ∂g dg .

(A1)

Knowing that ∂C/∂k x = -γC/(2k) and ∂C/∂g = C/(2g ), then

from ( 15). The x derivatives are then expanded and expressed as a function derivatives of C only. Especially, from (A1),

and

as dk x = γdk. After collecting the ∂ x C terms in (A2) using (A3) and (A4), eq. ( 17) is obtained.

The initial condition on C can be derived from the one on k x , equation (23):

where

and

A2 Equation for N Although equation ( 14) is already expressed as a function of N and C, it can be expressed in a more convenient way by expanding the x-derivative term. Especially:

Then, using (A4) and collecting terms in ∂ x C, equation ( 16) is obtained.

Appendix B Numerical Scheme

B1 General Procedure Any of the above equations can be translated in the frame moving with the long wave phase speed through a change of coordinates from (x, t) to (ψ, T ) where ψ is the long wave phase and T = σ L t accompanied with a change of variables from N and k x to q = (α, κ x ), where α = N (ψ, T )/N i (B1)

Equations ( 14) and ( 15) are then written in a nonlinear advection form

The flux function f is split into a quasi-linear term g and a ψ-dependent source term s f (q, ψ) = g(q) + s(q, ψ)

and

which includes the contributions of the orbital currents and the apparent gravity. Here, both functions are nonlinear, which requires a careful choice of numerical methods. We verified numerically that the scheme is conservative, in the sense that both the total action and wave number over the whole numerical domain (a long wave phase) are conserved during the time integration, within the round-off error.

The combination of a Runge-Kutta order 4 and the MUSCL4 numerical scheme of [START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF] is found to provide satisfying results for the numerical solution of the advection-only system and for the stability of the full solution computation. The Runge-Kutta scheme is used to solve the ODE

where discretized quantities are indexed with j and F n j±1/2 corresponds to MUSCL numerical fluxes. The latter are computed according to

Appendix C Multiple Time Scales Modulated Fields at Order ε L

Here we look for spatially periodic solutions for the short waves field by switching from the set of equations ( 14) and ( 15) on N and k x to the alternative set ( 14) and ( 17) on N and C, the equation on C being more convenient algebraically than the one on k x . Initial conditions ( 22) and ( 24) must be associated with the present set of equations. The variables are then expanded in powers of the small parameter ε L ,

Similarly for γ, which is a function of C only, we have:

At order ε 0 L , equations ( 16) and ( 17) and their initial conditions are

which are nonlinear forced advection equations. For simplicity, we consider the ε 0 L solutions to be constant throughout the propagation, although this system of equations would allow for more complicated solutions. Hence we consider:

and consequently

At order ε L , the equations to be solved are

because U is a ε L order function, according to (8).

In the present case, the solution of the system (C11) -(C14) would invalidate the perturbative expansion performed so far. Indeed, as it is discussed in section 33.2, this system presents a secular term -N 0 (3γ 2 0 -2)/(2γ 0 )∂ x C 1 in equation (C11). Its analytical solution would result in an unbounded growth of the solution N 1 namely, or instability. As a consequence, at large enough times, N 1 would become greater than N 0 , which would be in contradiction with the assumptions behind the perturbative expansion (C1). This was not considered in LH&S.

As an alternative to the perturbative expansion performed until now, and in the presence of a secular term, a multiple scale method can be employed (?, ?, e.g.)]Nayfeh1973 to revisit the problem from the beginning. This methods consists in keeping the perturbative expansion (C1) and (C2), and completing it by replacing the single time dimension by two. Two time scales are introduced, a fast one, t 0 = t and a slow one, t 1 = ε L t. t 0 is the scale of the regular oscillations, while t 1 is the scale of the instability. In the multiple scale method, these two variables are considered as independent. This introduces a single modification in the previous derivations. The time derivative is changed into

At order ε 0 L , equations (C4) to (C7) remain true only by replacing t by t 0 . The solutions (C8) and (C9) are kept, but the dependency on t 1 is not solved. In other words

At order ε 1 L , changes appear due to the introduction of the instability time. t 1 derivatives of N 0 and C 0 respectively are added in equations ( C11) and (C12), according to (C15). In the multiple time scales method, the slow time derivative terms are here to balance the secular terms, so that we are left with, on the one side,

and on the other,

with still the same initial conditions (C13) and (C14). The solutions are obtained by solving for (C18) and (C19) alone first, and introducing their solutions into the instability equations ( C21) and (C22).

The full solutions are summarized here. The phase speed consists in an average constant value

modulated with a fixed amplitude

where the constants χ and µ are defined by

and where the initial modulation on C is given by (A9). The action is already distorted at order ε 0 L :

where

with an additional modulation at order ε 1 L :


Coline Monchanin
email: coline.monchanin@univ-tlse3.fr

Erwann Drujont

Jean-Marc Devaud

Mathieu Lihoreau

Andrew B Barron

Metal pollutants have additive negative effects on honey bee cognition

Keywords: Apis mellifera, PER conditioning, pollutant interaction, arsenic, lead, copper

Honey bees displayed reduced learning and memory performance following acute exposure to arsenic, copper or lead. Exposure to combinations of these metals had additive effects.

Introduction

Metal pollution is of increasing concern for both ecosystem and public health [START_REF] Nriagu | Quantitative assessment of worldwide contamination of air, water and soils by trace metals[END_REF]. Over the last century, the widespread use of metals in domestic, industrial and agricultural applications [START_REF] Bradl | Sources and origins of heavy metals[END_REF] has considerably elevated their concentrations in water [START_REF] Mance | Pollution threat of heavy metals in aquatic environments[END_REF] and terrestrial habitats [START_REF] Krämer | Metal hyperaccumulation in plants[END_REF][START_REF] Su | A review on heavy metal contamination in the soil worldwide: situation, impact and remediation techniques[END_REF] to potentially toxic levels.

Pollinators, such as honey bees, are directly exposed to metal pollutants when foraging on contaminated nectar and pollen [START_REF] Perugini | Heavy metal (Hg, Cr, Cd, and Pb) contamination in urban areas and wildlife reserves: honeybees as bioindicators[END_REF][START_REF] Xun | Heavy metals in nectar modify behaviors of pollinators and nectar robbers: consequences for plant fitness[END_REF], and while flying through air containing suspended particles [START_REF] Thimmegowda | A field-based quantitative analysis of sublethal effects of air pollution on pollinators[END_REF]. Metals accumulate in the bodies of adults [START_REF] Giglio | Apis mellifera ligustica, Spinola 1806 as bioindicator for detecting environmental contamination: a preliminary study of heavy metal pollution in Trieste, Italy[END_REF] and larvae [START_REF] Balestra | Bees, honey, larvae and pollen in biomonitoring of atmospheric pollution[END_REF], as well as in the hive products [START_REF] Satta | Combination of beehive matrices analysis and ant biodiversity to study heavy metal pollution impact in a post-mining area (Sardinia, Italy)[END_REF]. For instance, concomitant bioaccumulation of arsenic (As), copper (Cu) and lead (Pb), resulting from metal production industries [START_REF] Kabir | Current status of trace metal pollution in soils affected by industrial activities[END_REF] and mining [START_REF] Khaska | Arsenic and metallic trace elements cycling in the surface water-groundwater-soil continuum down-gradient from a reclaimed mine area: isotopic imprints[END_REF][START_REF] Lee | Human risk assessment of As, Cd, Cu and Zn in the abandoned metal mine site[END_REF], is common in honey bees [START_REF] Badiou-Bénéteau | Honeybee biomarkers as promising tools to monitor environmental quality[END_REF][START_REF] Giglio | Apis mellifera ligustica, Spinola 1806 as bioindicator for detecting environmental contamination: a preliminary study of heavy metal pollution in Trieste, Italy[END_REF][START_REF] Goretti | Heavy metal bioaccumulation in honey bee matrix, an indicator to assess the contamination level in terrestrial environments[END_REF] and their honey [START_REF] Pisani | Minor and trace elements in different honey types produced in Siena County (Italy)[END_REF][START_REF] Terrab | Contribution to the study of avocado honeys by their mineral contents using inductively coupled plasma optical emission spectrometry[END_REF].

The deleterious effects of metals on humans [START_REF] Tchounwou | Heavy metal toxicity and the environment[END_REF] and some model animals (mice: Cobbina et al., 2015;flies: Do anlar et al., 2014) are well-known. As, Cu, Pb and other metals have neurotoxic effects that induce neural and neuromuscular alterations, sensory impairments and many related forms of behavioural dysfunctions [START_REF] Chen | Metals and neurodegeneration[END_REF]. Deficits in cognition and memory have been reported for As (e.g. humans: [START_REF] Tolins | The developmental neurotoxicity of arsenic: cognitive and behavioral consequences of early life exposure[END_REF]mice: Tyler et al., 2018;[START_REF] Wu | Effects of high fluoride and arsenic on brain biochemical indexes and learning-memory in rats[END_REF], Pb (e.g. mice: [START_REF] Anderson | Effects of low level lead exposure on associative learning and memory in the rat: Influences of sex and developmental timing of exposure[END_REF]humans: Mason et al., 2014) and Cu (e.g. mice: [START_REF] Lamtai | Chronic copper exposure leads to hippocampus oxidative stress and impaired learning and memory in male and female rats[END_REF][START_REF] Pal | Biochemical, histological, and memory impairment effects of chronic copper toxicity: a model for non-wilsonian brain copper toxicosis in wistar rat[END_REF]flies: Zamberlan, 2020).

Recent studies showed that low doses of Pb (Monchanin et al., 2021a) and selenium (Se) [START_REF] Burden | Acute exposure to selenium disrupts associative conditioning and long-term memory recall in honey bees (Apis mellifera)[END_REF] also impair behaviour and cognition in honey bees, suggesting a widespread impact on pollinators. So far, however, very little attention has been given to the potential combined effects of co-exposure to different metals [START_REF] Monchanin | Current permissible levels of heavy metal pollutants harm terrestrial invertebrates[END_REF].

Interactions among stressors are commonly classified as antagonistic (when the effect of one stressor reduces the effect of the other one), additive (when stressors have simple cumulative effects) or synergistic (when stressors together have a greater effect than the sum of their individual effects) [START_REF] Folt | Synergism and antagonism among multiple stressors[END_REF]. Additive effects of As, Cu and Pb have been described for humans [START_REF] Lin | Combined toxicity of heavy metal mixtures in liver cells[END_REF], rats [START_REF] Aktar | Individual and combined effects of arsenic and lead on behavioral and biochemical changes in mice[END_REF][START_REF] Mahaffey | Concurrent exposure to lead, cadmium, and arsenic: effects on toxicity and tissue metal concentrations in the rat[END_REF][START_REF] Schmolke | Arsenic-copper interaction in the kidney of the rat[END_REF] and fishes [START_REF] Verriopoulos | Combined toxicity of copper, cadmium, zinc, lead, nickel, and chrome to the copepod Tisbe holothuriae[END_REF]. In rats, for instance, co-exposure to Pb and As disrupted brain biogenic amine levels [START_REF] Agrawal | Changes in tissue oxidative stress, brain biogenic amines and acetylcholinesterase following co-exposure to lead, arsenic and mercury in rats[END_REF]. In humans, it was hypothesized that combined exposure to Pb and As, or other metal pollutants, have additive or synergistic toxic responses leading to cognitive dysfunction [START_REF] Karri | Heavy metals (Pb, Cd, As and MeHg) as risk factors for cognitive dysfunction: a general review of metal mixture mechanism in brain[END_REF]. To our knowledge, two studies have addressed the impact of metallic cocktails on bee physiology.

Honey bees simultaneously exposed to Pb, cadmium (Cd) and Cu accumulated significant levels of these metals in their body and had lower brain concentrations of dopamine compared to unexposed honey bees [START_REF] Nisbet | Preventive action of zinc against heavy metals toxicity in honeybee[END_REF]. Cd and Cu exerted a weak synergistic effect on honey bee survival [START_REF] Di | Joint effects of cadmium and copper on Apis mellifera foragers and larvae[END_REF]. However, none of these studies investigated potential effects of combined exposure on cognition.

Here we compared the effects of an exposure to single metals or ecologically relevant combinations of these metals on honey bee learning and memory. We hypothesised that combinations of metals may have synergistic negative effects, as it has been found with pesticides [START_REF] Yao | Influences of acephate and mixtures with other commonly used pesticides on honey bee (Apis mellifera) survival and detoxification enzyme activities[END_REF]Zhu et al., 2017). We tested individual honey bees in a standard protocol of proboscis extension reflex (PER) conditioning following acute exposure to As, Pb and Cu or a combination of them. We tested three concentrations of As, considered the most toxic substance (ATSDR, 2019), and added one concentration of Cu or Pb (binary mixtures), or both (tertiary mixture), to reach the molarity of the As solutions, allowing us to better assess any combined effects.

Materials and methods

Metal solutions

Arsenic (NaAsO2), lead (PbCl2) and copper (CuCl22H2O) were purchased from Sigma-Aldrich Ltd (Lyon, France) and diluted in 50% (w/v) sucrose solution. Control honey bees were fed 50% sucrose solution. Three concentrations of As were used (Table 1): a low concentration (0.13 µM) corresponding to the maximal permissible value in drinking water (0.01 mg.L -1 ) (Codex Alimentarius, 2015), a high concentration (0.67 µM) corresponding to half the maximal permissible value in irrigation water (0.1 mg.L -1 ) [START_REF] Ayers | Water quality for agriculture. Food and Agriculture Organization of the United Nations[END_REF], and an intermediate concentration (0.40 µM). This range of concentrations was reported in water sampled from polluted areas (e.g. mining sites) and in honey (Table S1). For Pb and Cu, we chose 0.27 µM (0.055 mg.L -1 of Pb and 0.017 mg.L -1 of Cu) so that the binary combinations (As 0.13 µM + Cu 0.27 µM or As 0.13 µM + Pb 0.27 µM) could be compared to the As intermediate concentration (0.40 µM), and the tertiary combination (As 0.13 µM + Pb 0.27 µM + Cu 0.27 µM) to the As high concentration (0.67 µM) (Table 1). These concentrations of Pb and Cu have also been reported in honey samples (Table S1). The mass consumed for As and the concentrations for Cu and Pb fell within sublethal ranges for the honey bee: the LD50 of elemental As for NaAsO2 ranged from 0.330 to 0.540 µg/bee [START_REF] Fujii | Oral dose toxicity vs. tissue residue levels of arsenic in the honey bee (Apis mellifera L.)[END_REF], the LC50 of Cu is 72 mg.L -1 [START_REF] Di | Laboratory bioassays on the impact of cadmium, copper and lead on the development and survival of honeybee (Apis mellifera L.) larvae and foragers[END_REF] and of Pb is 345 mg.L -1 [START_REF] Di | Laboratory bioassays on the impact of cadmium, copper and lead on the development and survival of honeybee (Apis mellifera L.) larvae and foragers[END_REF].

Bee exposure to metals

We collected honey bees (Apis mellifera) returning from foraging trips at the entrance of five different hives in mornings during August 2020. We anesthetised the bees on ice and harnessed them in plastic tubes, secured with tape and a droplet of wax at the back of the head [START_REF] Matsumoto | Revisiting olfactory classical conditioning of the proboscis extension response in honey bees: A step toward standardized procedures[END_REF]. We tested all the bees for intact proboscis extension (PER) by stimulating their antennae with 50% sucrose. We then fed the responding honey bees 5 µL of 50% sucrose solution (see Table 1), making sure they consumed the whole droplet, and left them to rest for 3 h in the incubator (temperature: 25±2°C, humidity: 60%). Honey bees that did not respond to the sucrose solution were discarded.

Absolute learning

Prior to conditioning, we tested all honey bees for PER by stimulating their antennae with 50% sucrose solution, and kept only those that displayed the reflex. We then performed olfactory absolute conditioning according to a standard protocol using an automatic stimulus delivery system [START_REF] Aguiar | Can honey bees discriminate between floral-fragrance isomers?[END_REF]. Honey bees had to learn to respond to an olfactory conditioned stimulus (CS, 1-nonanol, Sigma-Aldrich Ltd, Lyon, France) reinforced with unconditioned stimulus (US, 50% sucrose solution), over five conditioning trials with a ten-minute inter-trial interval. Each trial (37 s in total) began when a bee was placed in front of the stimulus delivery system, which released a continuous flow of clean air (3,300 mL.min -1 ) to the antennae. After 15 s, the odour was introduced into the airflow for 4 s, the last second of which overlapped with sucrose presentation to the antennae using a toothpick. This was immediately followed by feeding for 4 s by presenting the toothpick to the proboscis. The bee remained for another 15 s under the clean airflow. We recorded the presence or absence (1/0) of a conditioned PER in response to the odorant presentation during each conditioning trial. Honey bees spontaneously responding in the first conditioning trial were discarded from the analysis. The sum of conditioned responses over all trials provided an individual acquisition score (between 0 and 4), and honey bees responding at the last trial were categorized as learners.

Long-term memory

Only honey bees that had learnt the task were kept for the analysis of memory performance.

After conditioning, these honey bees were fed 15 µL of 50% sucrose solution, left overnight in the incubator, and fed another 5 µL of sucrose solution the following morning. Three hours after (24 h post-conditioning), we performed the retention test, consisting of three trials similar to conditioning except that no sucrose reward was presented. In addition to the odour used during the conditioning (CS), we presented two novel odours, in randomized order, to assess the specificity of the memory: nonanal was expected to be perceived by honey bees similarly to 1-nonanol, while 1-hexanol was expected to be perceived differently [START_REF] Guerrieri | Perceptual and neural olfactory similarity in honeybees[END_REF].

We recorded the presence or absence (1/0) of a conditioned PER to each odorant at each memory retention trial. We classified honey bees according to their response patterns: response to the CS only, response to the CS and the similar odour (low generalization level), response to all odours (high generalization level), no or inconsistent response.

Statistics

We analysed the data using R Studio v.1.2.5033 (RStudio Team, 2015). Raw data are available in Dataset S1. We performed binomial generalised linear mixed-effects models (GLMM) (package lme4; [START_REF] Bates | Fitting linear mixed-effects models using lme4[END_REF], with hive and conditioning date as random factors and treatment as fixed effect. Using the GLMMs, we evaluated whether molarity or treatment impacted the initial response to antennal stimulation, the spontaneous response in the first conditioning trial, the response in the last trial, the response to each odorant during the memory test, the proportion of honey bees per response pattern in the retention test, and the survival at 24 h. Acquisition scores were standardised and compared with GLMMs using Template Model

Builder [START_REF] Brooks | glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling[END_REF]. For all response variables, we compared (1) the treated groups to the control, (2) groups exposed to concentrations of the same molarity (e.g. Med [As], [As+Cu] and [As+Pb]), (3) the separate and joint effects of the treatments (e.g. Low [As], [Cu] and

[As+Cu]) in order to identify interactive effects (antagonistic, additive, synergistic).

Results and discussion

Exposure to metals did not impact appetitive motivation

The proportion of honey bees that responded to the initial antennal stimulation with sucrose was similar among treatments (GLMM: p>0.05). All honey bees responded in the other groups.

Therefore, treatment did not affect appetitive motivation or sucrose perception. Consistently with our observations, the ingestion of similar concentrations of Pb and Cu had no effect on the responsiveness to increasing concentrations of sucrose [START_REF] Burden | Acute sublethal exposure to toxic heavy metals alters honey bee (Apis mellifera) feeding behavior[END_REF]. By contrast, [START_REF] Di | Joint effects of cadmium and copper on Apis mellifera foragers and larvae[END_REF] found that honey bees exposed to increasing concentrations of a mixture of Cu and Cd exhibited a decreased ability to distinguish sucrose concentrations, but this may be explained by the much higher (at least 600 times) concentrations used in that study. Thus, in our conditions any impact on appetitive learning is unlikely due to a decreased motivation for sucrose.

Individual and joint exposures to metals reduced learning performance

Two out of the 381 honey bees submitted to the absolute learning task spontaneously responded to the first odour presentation and were therefore discarded. In all groups, the number of honey bees showing the conditioned response increased over trials, thus showing learning (Fig. 1A).

However, fewer honey bees exposed to metals learned the task when compared to controls (GLMM: p<0.05, except for Low [As], p=0.082). Accordingly, the acquisition scores of honey bees from all treatments were lower than controls (Fig. 1B). Honey bees exposed to Med

Individual and joint exposures to metals reduced long-term memory specificity

To examine possible effects of metal exposure on memory retention, we tested the capacity for long-term memory retention. Only honey bees that had learned the CS-US association at the end of conditioning were tested, to evaluate retention levels independently of possible biases due to slight variations in learning performance among treatments. 167 out of the 379 honey bees submitted to the absolute learning task did not learn and were therefore discarded. We found no effect of treatment on survival at 24 h (GLMM: p>0.05). However, long-term memory was significantly affected (Fig. 2). Overall, treated honey bees responded less to the learned odorant (CS) than controls, as indicated by a significant effect of exposure to metals on retention levels (GLMM: p<0.05) (Fig. 2A). Yet, this decrease was not significant for honey bees exposed to Med [As] (GLMM: -0.260±0.628, p=0.679) and High [As] (GLMM: -1.023±0.570, p=0.073). Finally, there was no clear dose effect on responses to the CS among treated groups (GLMM: -0.576±0.579, p=0.320).

Yet, individual response patterns (Fig. 2B) revealed a loss of memory specificity. While honey bees from all treatments responded equivalently to the similar odour (GLMM: p>0.05), those exposed to higher doses responded more frequently to all odorants, indicating a higher degree of response generalization (GLMM: 1.954±0.775, p=0.012). This was accompanied by a significantly lower proportion of specific (CS-only) responses for honey bees exposed to [Pb] (GLMM: -1.795±0.690, p=0.009), low [As] (GLMM: -1.313±0.589, p=0.026) and

[As+Cu+Pb] (GLMM: -1.200±0.588, p=0.041). Exposure also significantly increased the frequency of inconsistent responses as compared to controls (GLMM: p<0.05). This was the case for each individual treatment except for Med [As], p=0.293). Thus, exposure to metals had a negative impact on memory performance at 24h. The analysis of individual response patterns also revealed additive effects as they neither differ among groups exposed to solutions with the same molarity, nor between single and mixed metal treatments (GLMM: p>0.05). Thus, most treatments reduced memory performance at 24h.

The additive effects of metal mixtures may be explained by common pathways of toxicity

Although many mechanisms of metal toxicity have not yet been elucidated, some points of consensus are emerging from the literature. Firstly, interactions between metals can occur outside the organism (including interactions with the environment [START_REF] Grobelak | Heavy metal mobility in soil under futuristic climatic conditions[END_REF][START_REF] Noyes | Forecasting the impacts of chemical pollution and climate change interactions on the health of wildlife[END_REF]), and during uptake into the organism, leading to potentially toxic processes of speciation, absorption, binding, transport and distribution [START_REF] Wu | A review of toxicity and mechanisms of individual and mixtures of heavy metals in the environment[END_REF].

Once metals enter an organism, they can induce, interact or inhibit a range of biological responses and metabolic pathways. By mimicking other essential metals [START_REF] Bridges | Molecular and ionic mimicry and the transport of toxic metals[END_REF] or damaging the permeability of biological membranes [START_REF] Rothshein | Cell membrane as site of action of heavy metals[END_REF], metals allow the uptake or loss of other compounds in/out intracellular compartments [START_REF] Viarengo | Heavy metal cytotoxicity in marine organisms: effects on Ca2+ homeostasis and possible alteration of signal transduction pathways[END_REF].

Metals are known to disrupt signalling and calcium homeostasis (particularly important in neurons) by interfering with the calcium channels [START_REF] Bridges | Molecular and ionic mimicry and the transport of toxic metals[END_REF][START_REF] Chavez-Crooker | Copper transport by lobster hepatopancreatic epithelial cells separated by centrifugal elutriation: measurements with the fluorescent dye Phen Green[END_REF][START_REF] Tamano | Dynamic action of neurometals at the synapse[END_REF]. This might lead to dysfunction and cytotoxicity due to the disruption of cell signalling and calcium homeostasis. Genotoxicity (Do anlar et al., 2014) may be achieved through covalent binding to DNA [START_REF] Brocato | Basic mechanics of DNA methylation and the unique landscape of the DNA methylome in metal-induced carcinogenesis[END_REF][START_REF] Senut | Lead exposure disrupts global DNA methylation in human embryonic stem cells and alters their neuronal differentiation[END_REF]. Eventually, oxidative stress and lipid peroxidation of the cell membrane may lead to neuronal death. Additionally, metals in mixture could interact at target sites, but the potential for modulation of that toxic impact is largely unknown [START_REF] Svendsen | Toxicokinetics and toxycodynamics[END_REF]. Metal mixtures could change the bioavailibity, toxicokinetic and toxicodynamic of each metal, that could directly impact the toxicity towards the organism [START_REF] Løkke | Tools and perspectives for assessing chemical mixtures and multiple stressors[END_REF]. Based on these shared mechanisms of toxicity that include oxidative stress [START_REF] Nikoli | The impact of sublethal concentrations of Cu, Pb and Cd on honey bee redox status, superoxide dismutase and catalase in laboratory conditions[END_REF][START_REF] Zaman | An insect model for assessing oxidative stress related to arsenic toxicity: arsenic toxicity and oxidative stress[END_REF], apoptosis [START_REF] Raes | Copper induces apoptosis in Aedes C6/36 cells[END_REF] and interference with neurotransmitters [START_REF] Nisbet | Preventive action of zinc against heavy metals toxicity in honeybee[END_REF], the toxic effects of metal pollutants in mixtures may be expected to be additive [START_REF] Von Stackelberg | Mixtures, metals, genes and pathways: a systematic review[END_REF].

These effect may affect many aspects of neural activity and brain function in honey bees, as in other species [START_REF] Karri | Heavy metals (Pb, Cd, As and MeHg) as risk factors for cognitive dysfunction: a general review of metal mixture mechanism in brain[END_REF]. Here, we focused on learning and memory of olfactory cues because they play crucial roles in the behavioural ecology of honey bees and other pollinators, for the identification of food resources. Our results in controlled laboratory conditions suggest that exposure to sublethal combinations of toxic elements in the field might alter individual foraging efficiency, and in turn jeopardize survival of pollinator populations.

Our findings call for further evaluation of the joint actions of metals [START_REF] Meyer | Metal mixtures modeling evaluation project: background[END_REF] to better assess the risk they pose [START_REF] Nys | A framework for ecological risk assessment of metal mixtures in aquatic systems: Ecological risk assessment framework for metal mixtures[END_REF][START_REF] Otitoloju | Relevance of joint action toxicity evaluations in setting realistic environmental safe limits of heavy metals[END_REF] and better inform regulatory framework (European Commission, 2012). Current risk assessment guidance mainly assesses the effect of individual exposure, that fails to capture potential interactive effects. Hence, the evaluation of metal mixture impacts and their modes of action needs to be developed [START_REF] Sasso | A generalized physiologicallybased toxicokinetic modeling system for chemical mixtures containing metals[END_REF]. More generally, the study of the interactions between toxic metals and environmental factors [START_REF] Naqash | Interaction of freshwater microplastics with biota and heavy metals: a review[END_REF] and their impacts on the toxicokinetic of other chemicals (EFSA Scientific Committee et al., 2019) (e.g. pesticides [START_REF] Sgolastra | Lethal effects of Cr(III) alone and in combination with propiconazole and clothianidin in honey bees[END_REF][START_REF] Singh | Synergistic effects of heavy metals and pesticides in living systems[END_REF], volatile organic compounds [START_REF] Sasso | A generalized physiologicallybased toxicokinetic modeling system for chemical mixtures containing metals[END_REF] etc.) should be implemented in this integrated research framework.

Conclusion

In summary, we demonstrated that arsenic, lead, copper or combinations of these metals, at levels found in the environment, slow down appetitive learning and reduce long-term memory specificity in honey bees. These metals show simple additive effects as we found no differences between different solutions of the same molarity suggesting possible non-linear effects (synergism or antagonism). Thus, regarding effects on learning and memory, concentration seem to be more important than identity of any specific metal. Since learning and memory of olfactory cues play crucial roles in the behavioural ecology of honey bees, acute exposure to metal pollutants mixtures could impair fundamental hive function and population growth.

  [As] (GLMM: -0.610±0.246, p=0.013), High [As] (GLMM: -0.639±0.241, p=0.008) and[As+Cu+Pb] (GLMM: -0.592±0.244, p=0.015) had acquisition scores significantly lower than controls. Honey bees exposed to solutions of [As+Pb] had similar acquisition scores to Med[As] (GLMM: 0.299±0.234, p=0.201), but honey bees exposed to [As+Cu] performed better (GLMM: 0.596±0.241, p=0.013). Honey bees exposed to High [As] and [As+Cu+Pb] exhibited similar acquisition scores (GLMM: p=0.810). We found no difference in the acquisition scores and the proportions of learners between honey bees treated with a single metal and mixed treatments (GLMM: p>0.05), that would have indicated non-additive effects (i.e. antagonistic or synergistic). Thus, exposure to metals significantly reduced learning performance, and combined exposure appeared to exert simple additive deleterious effects.

Figure 1 :

 1 Figure 1: Learning. A) Learning curves show changes in the percentages of honey bees displaying the conditioned proboscis extension response (PER) over five training trials. Asterisks indicate significant differences in responses at the last trial compared to control honey bees. B) Violin plots of acquisition score values (sum of conditioned responses for each honey bee). Symbols (circle: single exposure; triangle: binary mixture; diamond: tertiary mixture) indicate the mean score for each treatment. Significant differences between groups exposed to the same molarity solutions (#) or with respect to control honey bees (*) are indicated (#/*p<0.05, **p<0.01; GLMM).

Figure 2 :

 2 Figure 2: Long-term memory. A) Percentages of responses to the CS odour in the 24 hmemory retention test (mean ± s.e.m). B) Distribution of honey bees according to their individual response pattern during the long-term memory test: response to CS only; response to CS and similar; response to all odours; no or inconsistent response. Significant differences with controls are indicated (*p<0.05, **p<0.01; GLMM).

Table 1 :

 1 Concentrations used. Combined treatments are shown in grey.

			As	Cu	Pb	As	Cu	Pb
	Control	0	0	0	0	0	0	0
	Low [As]	0.13	0.01	0	0	0.05	0	0
	[Cu]	0.27	0	0.02	0	0	0.09	0
	[Pb ]	0.27	0	0	0.06	0	0	0.28
	Med [As]	0.40	0.03	0	0	0.15	0	0
	[As+Cu]	0.40	0.01	0.02	0	0.05	0.09	0
	[As+Pb]	0.40	0.01	0	0.06	0.05	0	0.28
	High [As]	0.67	0.05	0	0	0.25	0	0
	[As+Cu+Pb] 0.67	0.01	0.02	0.06	0.05	0.09	0.28

Treatment

Molarity (µM) Concentration (mg.L -1 ) Ingestion of 5µL (ng/bee)

Acknowledgments

We thank Olivier Fernandez for assistance with beekeeping.

Funding

This work was supported by the CNRS. CM was funded by a PhD fellowship from French Ministry of Higher Education, Research and Innovation. ABB was funded by a Future Fellowship from the Australian Research Council (FT140100452) and the Eldon and Anne Foote Trust. ML was funded by grants of the Agence Nationale de la Recherche (ANR-16-CE02-0002-01, ANR-19-CE37-0024, ANR-20-ERC8-0004-01) and the European Regional Development Fund (project ECONECT).

Data availability

Raw data will be available on Dryad repository upon publication.

Competing interests

The authors declare no competing or financial interests.

Zamberlan, D. C. (2020).

Copper decreases associative learning and memory in Drosophila melanogaster. Sci. Total Environ. 9, 1-14. Zhu, Y. C., Yao, J., Adamczyk, J. and Luttrell, R. (2017). Synergistic toxicity and physiological impact of imidacloprid alone and binary mixtures with seven representative pesticides on honey bee (Apis mellifera). PLOS ONE 12, e0176837.