Mathieu Chatelain
email: chatelain@obs-banyuls.fr

Katell

Katell Guizien

Modelling coupled turbulence -dissolved oxygen dynamics near the sediment-water interface under wind waves and sea swell

Keywords:

come

Modelling coupled turbulence -Dissolved oxygen dynamics near the sediment-water interface under wind waves and sea swell

Introduction

Diffusional mass transfer between bottom sediments and overlying water constitutes an essential coupling for benthic and pelagic ecosystems [START_REF] Boudreau | The benthic boundary layer: Transport processes and biogeochemistry[END_REF]. Transfer of oxygen across the sediment-water interface (SWI) regulates biological and geochemical processes of organic matter degradation in the upper sediment. Vertical oxygen mass transfers at the benthic or bottom boundary layer occur in a zone where gradients in physical, chemical and biological properties are sharp. In the water column, solute transport is dominated by turbulent mixing, except within the diffusive boundary layer (DBL) where turbulence diffusivity becomes negligible compared to oxygen molecular diffusion [START_REF] Gundersen | Microstructure of diffusive boundary layers and the oxygen uptake of the sea floor[END_REF]. The thickness of the DBL regulates the kinetics of oxygen supply to the benthic organisms: while a thick DBL associated with small steady flow velocities (laminar regime) may be limiting for organisms with high oxygen uptake rates, a thin DBL associated with large steady flow velocities (turbulent regime) results in the highest possible oxygenation [START_REF] Nakamura | Effect of flow velocity on sediment oxygen demand: Theory[END_REF][START_REF] Hondzo | Dissolved oxygen transfer at the sediment-water interface in a turbulent flow[END_REF][START_REF] Steinberger | Diffusional mass transfer at sediment-water interface[END_REF]. In nature, flows are generally turbulent with DBL thicknesses on the order 1-2 mm [START_REF] Jørgensen | Diffusive boundary layers and the oxygen uptake of sediments and detritus[END_REF] but they are also more often unsteady. As a consequence, the DBL thickness will fluctuate inducing oxygen profile dynamics [START_REF] Jørgensen | The diffusive boundary layer of sediments: Oxygen microgradients over a microbial mat[END_REF][START_REF] Glud | Effect of the diffusive boundary layer on benthic mineralization and O2 distribution: A theoretical model analysis[END_REF]. Yet, such dynamics are not instantaneously responding to the hydrodynamic forcing. Transient oxygen profiles were observed in sediment cores immediately after stirring was stopped [START_REF] Lohse | Oxygen pore water profiles in continental shelf sediments of the North Sea: Turbulent versus molecular diffusion[END_REF]. Numerical simulations have shown that the duration of the transitory regime between two steady states (with and without flow) decreases as both oxygen consumption and flow velocity increase [START_REF] Higashino | Unsteady diffusional mass transfer at the sediment/water interface: Theory and significance for SOD measurement[END_REF].

Amongst unsteady flows, oscillatory flows are of particular interest for the diffusion-consumption of oxygen at the SWI. Periodic flows are frequent in nature as they are linked to surface gravity waves.

In the oceans, such flows are the most energetic [START_REF] Massel | Ocean surface waves. Their physics and prediction[END_REF] and they often destabilize the SWI in the nearshore zone [START_REF] Harris | Estimating global shelf sediment mobility due to swell waves[END_REF]. Surface gravity wave energy is equally distributed between two types characterized by their period: wind waves and swell (periods from 2 to 30 s) and the diurnal and semi-diurnal tides (periods from 12 to 24 h). In lakes and lagoons, the period range for surface gravity waves is generally restricted to wind waves (period < 8 s) due to fetch limitations. But other types of oscillatory flows like a seiche with even longer periods (depending on the basin lengthscale and stratification) can be observed [START_REF] Proudman | Dynamical Oceanography[END_REF]. In 2003, the oxygen profile dynamics linked to periodic bottom layer turbulence under low-frequency waves (lake internal seiching with period of 18 hours) was observed for the first time [START_REF] Lorke | Breathing sediments: The control of diffusive transport across the sediment-water interface by periodic boundary-layer turbulence[END_REF]. Additionally, the effect of periodic flows (periods ranging from minutes to hours) on oxygen diffusion in the water column was investigated by numerical simulations [START_REF] Higashino | Periodic diffusional mass transfer near sediment/water interface: Theory[END_REF].

The present study examines the effect of wind waves and sea swell (periods ranging from 2 to 30 s) on oxygen diffusion and consumption at the SWI. Oscillatory flows exhibiting large velocity fluctuations from zero (at flow reversal) to a maximum value (orbital velocity) are associated with these highfrequency gravity waves [START_REF] Airy | Tides and waves[END_REF]. Near the bed, an oscillatory bottom boundary layer develops which exhibits strong turbulence dynamics with flow relaminarisation around reversal [START_REF] Jensen | Turbulent oscillatory boundary layers at high Reynolds numbers[END_REF] and consequently, large DBL thickness fluctuations at the flow periodicity. Under oscillatory flows, the oxygen concentration distribution dynamics may result not only from the DBL thickness value, but also from its dynamics as well. Hence, a new one-dimensional vertical, unsteady numerical model for diffusion-consumption of dissolved oxygen (DO) above and below the SWI was driven by the turbulent diffusivity associated with the wave bottom boundary layer modelled by [START_REF] Guizien | 1DV bottom boundary layer modeling under combined wave and current: Turbulent separation and phase lag effects[END_REF]. Wind waves and sea swell were described by idealized monochromatic waves to assess the respective influence of the wave period and orbital velocity. 2 Material and methods

Model formulation

Dissolved oxygen traverses two diffusive boundary layers from the water column into the sediment.

The first layer is in the water immediately above the sediment, and the second is in the sediment immediately below the sediment-water interface (Fig. 1).

[Figure 1] A one-dimensional vertical (1DV) numerical model of the DO balance in water and sediment is proposed. For the sake of simplicity, we assume a constant porosity of the sediment. Advection due to permeability, bioturbation and bioirrigation in the sediment is neglected. Microbial organic matter degradation occurring in the sediment is taken into account, as well as chemical reactions involving oxygen consumption. Both are modelled by Michaelis-Menten consumption kinetics [START_REF] House | Factors influencing the extent and development of the oxic zone in sediments[END_REF]. Such a mathematical formulation allows a constant consumption when the DO concentration is large (i.e. non limiting, zero-order kinetics) and a decreasing consumption when DO concentration approaches zero (first-order kinetics). This formulation is thus consistent with both zero flux and concentration for DO at the lower boundary condition deep in the sediment. Moreover, prescribing a constant consumption near the interface assumes that neither acclimatization nor growth of the respiring community occurs inside a wave cycle. Finally, it is assumed that oxygen consumption can be neglected in the water column. Thus we define:

∂C * ∂t * = ∂ ∂z * D * ∂C * ∂z * -α. µ * C * K * O 2 + C * (1)
where the non-dimensional parameters are

C * = C C 0 , t * = νt h 2 , z * = z h , D * = D ν , µ * = µh 2 νC 0 , K O2 * = K O2 C 0 . ( 2
)
C is the DO concentration, t is time, z is the vertical coordinate (positive upward), D is the vertical diffusivity, µ is the maximum oxidation rate, K O2 is the oxygen half-saturation constant, and α is a numerical constant (α=1 in sediment and α=0 in water). Reference values are ν the kinematic viscosity in water, C 0 the bulk water concentration, and h a length greater than the sum of DBL thickness in water and DO penetration depth in sediment.

In sediment, oxygen vertical diffusivity D is equal to the effective diffusivity D s , which is molecular diffusion corrected for tortuosity. It can be expressed, through Archie's law, as a function of molecular diffusivity and porosity of the sediments [START_REF] Ullman | Diffusion coefficients in nearshore marine sediments[END_REF]:

D s = φ m-1 D w (3)
where φ is the porosity, and m denotes an exponent corresponding to different kinds of sediment (m=3, [START_REF] Manheim | Diffusimetry (diffusion constant estimation) on sediment cores by resistivity probe[END_REF]. Molecular diffusion D w for oxygen is assumed constant, using a Schmidt number Sc=ν/D w =500 at 20 o C [START_REF] Denny | Air and water: The biology and physics of life's media[END_REF].

In water, oxygen vertical diffusivity D is the sum of molecular (D w ) and turbulent (D t ) diffusivities.

Using the analogy between momentum and mass transfer, the turbulent diffusivity (D t ) is assumed to be equal to the eddy viscosity (ν t ) in water. The model formulation thus allows time-variations of turbulent diffusivity to account for unsteady hydrodynamics.

Equation ( 1) is solved using the implicit finite control volume method of [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. The spatial computational domain went from -h s in the sediment, to h w in the water column (h = h s + h w ).

Two exponential grids were defined with n s points in sediment and n w points in water. Each mesh starts at z=0, reading:

|z i+1 -z i | = dz i = s r i (4)
taking (r,s) equals (r s ,s s ) in sediment and (r w ,s w ) in water. The SWI is the first sediment point, located at z=0 (α=1, D = D s ). This grid was adopted because it is well adapted to refining the nearbed description of each domain where gradients are sharp. Meshing independence of the computations was studied and showed insensitivity of the results for a resolution that yields r s =1.054, s s =1×10 -5

(n s =194) and r w =1.054, s w =1×10 -7 (n w =320) for sediment and water, respectively (in Eq. 4).

Boundary conditions for DO concentration are imposed at the top and the bottom of the calculation domain:

C * = 0 at z * = -h s /h C * = 1 at z * = h w /h (5)
The first condition expresses that all DO reaching the sediment should be utilized within the oxygen penetration depth. The second indicates that at the upper limit of the computational domain the water column is fully mixed and the oxygen concentration equals the bulk water concentration. Bulk DO concentration was arbitrarily set to 10 g m -3 [START_REF] Glud | Seasonal dynamics of benthic oxygen uptake in a semienclosed bay: Importance of diffusion and faunal activity[END_REF].

As we are dealing with periodic hydrodynamic forcing, a finite number of time steps (n t =360) is required to define an entire period cycle, and periodic boundary conditions are applied stating C * (z,n t +1)=C * (z,1). The equilibrium solution is the solution when the maximum DO concentration difference at the same phase during two consecutive cycles was less than a convergence criterion . We used =10 -8 as the convergence criterion.

Sediment characteristics and hydrodynamic forcing

This study examines the nearbed dynamics of oxygen mass transfer transport under wave forcing over a smooth bottom, as can be observed in shallow estuaries open to the sea. Porous silty sediment conditions were used (median diameter d 50 =0.002 cm; bottom roughness k n =0.005 cm; porosity φ = 0.9). The half-saturation constant for DO is set to a constant value (K O 2 =0.2 g m -3 , [START_REF] Hao | The half-saturation coefficient for dissolved oxygen: A dynamic method for its determination and its effect on dual species competition[END_REF].

Two biogeochemical activities were tested: high oxidation rate (µ=2000 g m -3 d -1 ) and low oxidation rate (µ=50 g m -3 d -1 , [START_REF] Higashino | Unsteady diffusional mass transfer at the sediment/water interface: Theory and significance for SOD measurement[END_REF]. These values correspond to oxygen consumption times in the sediment of 422 s and 4.8 h, respectively.

A linear wave forcing is considered:

U (t) = U w sin 2πt T ( 6
)
where T is the leading wave period and U w is the nearbed orbital velocity.

In order to investigate the effect of oscillating flows typical of wind waves to swell on DO uptake dynamics, we considered wave orbital velocities ranging from 20 to 200 cm s -1 and periods ranging from 2 to 30 s. Virtually, any orbital velocity / period couple can be considered, provided the wave stability criterion is fullfilled. We excluded some orbital velocity / period couples which are unrealistic as they would correspond to unstable breaking waves. Using a breaking-wave criterion of H < 0.8 D (where H is the wave height and D is the water depth), the orbital velocity restrictions only applied to waves with periods shorter than 4 s.

Hydrodynamic forcing for the DO model consisted of time-varying eddy viscosity profiles ν t (z, t) in the wave boundary layer during a wave cycle. Together with bottom shear stress time series, these were computed using a 1DV bottom boundary layer numerical model dedicated to periodic flow [START_REF] Guizien | 1DV bottom boundary layer modeling under combined wave and current: Turbulent separation and phase lag effects[END_REF], given the bottom roughness k n . The shear velocity u * at the bed is then:

u * (t) = τ (t) ρ (7)
where ρ denotes the water density.

[Figure 2] Fig. 2 shows the time evolution during a wave cycle of the shear velocity u * (t) for two contrasted wave orbital velocities at the same wave period (T =15 s) to illustrate the laminar and turbulent regime.

In both regimes, during each half of a wave cycle, shear velocity varies from nearly zero to a maximum u * m (not a sine form). As typically seen in oscillatory flows, the shear velocity maximum and minimum are ahead of the wave velocity outside the boundary layer, showing that small velocities (near the bottom) reverse more easily than larger ones (far from the bottom) when the pressure gradient reverses. This phase shift decreases from 45 o in the laminar case to about 10 o in the fully turbulent regime [START_REF] Fredsøe | Mechanics of coastal sediment transport[END_REF]. For the same reason, turbulent diffusivity (not shown) exhibits a phase shift with the outer flow velocity which varies across the boundary layer. A laminar regime is defined as a wave period and an orbital velocity for which the phase shift is 45 o . It corresponds in our computations to a wave Reynolds number below 320:

Re δ = U w δ ν (8)
with δ = νT /π the wave viscous boundary layer height [START_REF] Vittori | Direct simulation of transition in an oscillatory boundary layer[END_REF]. In the turbulent regime, the shear velocity time evolution jumps during the acceleration phase, marking the onset of turbulence that occurs after the boundary layer relaminarisation around flow reversal [START_REF] Jensen | Turbulent oscillatory boundary layers at high Reynolds numbers[END_REF]. In the laminar regime, such a jump is absent (flow remains laminar during the entire wave cycle). The laminar regime (Re δ ≤ 320) is shown in light grey on Fig. 3. For the wave climates that were investigated, maximum shear velocities ranged from less than 1 to 8 cm s -1 . The sensitivity to the wave period is accentuated when the period is small. The dark grey area indicates breaking wave values on the same figure.

[Figure 3] Consequently, the oxygen DBL thickness δ D (t) (the area where D t < D w =2×10 -5 cm 2 s -1 ) also

shows temporal variations during a wave cycle, ranging from around 5×10 -3 to 3×10 -2 cm for T =25 s and U w =120 cm s -1 (u * m =5 cm s -1 ). It should be noted that the exponential mesh creates a coarser resolution in δ D (t) computations for large DBL thicknesses.

[Figure 4]

Model analysis

Computational results consisted of DO vertical profiles at different phases during a wave cycle.

For each DO profile, we defined the interfacial concentration C w (t), the penetration depth δ s (t) (with

C * (z=δ s (t))=0
.001) and the dissolved oxygen uptake (DOU) at the SWI (Fig. 1). The latter is evaluated from Fick's first law of diffusion [START_REF] Berner | Early diagenesis: A theoretical approach[END_REF]:

DOU (t) = -D w ∂C ∂z = - D w C 0 β(t) . ( 9
)
Practically, the instantaneous DOU(t) is computed taking the slope β(t) of the DO vertical profile over the three first grid points in water, which height is always much smaller than the DBL thickness. The dimensionless Sherwood number, Sh, is introduced, following [START_REF] Higashino | Periodic diffusional mass transfer near sediment/water interface: Theory[END_REF]:

Sh(t) = DOU (t) Sc u * m (C w (t) -C 0 ) (10)
where u * m is the maximum shear velocity. The Sherwood number can be interpreted as a ratio between the effective oxygen transfer rate and turbulent diffusion rate. It thus measures the turbulence efficiency to supply oxygen to the sediment. Following [START_REF] Lorke | Breathing sediments: The control of diffusive transport across the sediment-water interface by periodic boundary-layer turbulence[END_REF], we defined a time scale t D (t) for the oxygen diffusion across the DBL as:

t D (t) = δ D (t) 2 D w . ( 11
)
All these time-varying quantities are described by their mean value (indicated by the overline notation)

and their variance (VAR(x)=100×RMS(x)/x in %, where RMS is Root Mean Square) during a wave cycle.

To go further, a time scale (t w ) describing the temporal dynamics of the DBL thickness is introduced:

t w (t) = δ D (t) | ∂δ D ∂t | . ( 12
)
This time scale describing the intra-wave DBL dynamics should be compared to the time scale for the establishment of a local DBL. Similar reasoning was adopted by [START_REF] Lorke | Breathing sediments: The control of diffusive transport across the sediment-water interface by periodic boundary-layer turbulence[END_REF], where the authors compared relative importance of vertical diffusion versus horizontal advection. When t w (t) is larger than t D (t), the DBL thickness varies less rapidly than the time required for the DBL to establish: diffusive fluxes will follow the local DBL thickness fluctuations. Conversely when t w (t) is smaller than t D (t), the DBL thickness varies more rapidly than the time required for the DBL to establish: transient diffusive fluxes corresponding to average DBL thickness over the period are simulated.

Results

Intra-wave cycle dynamics

In the laminar regime (shear velocity time series shown on Fig.

2 for T =15 s, U w =10 cm s -1 , u * m =0.8 cm s -1 , Re δ =218
), the time required by oxygen to diffuse through the DBL (t D ) during the wave cycle was very large (≈ 1.1 days) compared to both the oxygen consumption time in the sediment (t c ) and the wave period (T ). The values of t D were also fairly constant during the wave cycle (data not shown). Consequently, both DOU and concentration at the SWI displayed very small and constant values (DOU =0.12 g m -2 d -1 , C w =0.01 C 0 ). The oxygen penetration depth was also very small (δ s =0.03 cm). When t D t c , oxygen is consumed more rapidly in the sediment than it diffuses from water into the sediment: SWI oxygenation is limited by the transfer time through the water DBL.

In the turbulent regime, the time required by oxygen to diffuse through the DBL (t D ) is much shorter than in the laminar regime, and reaches lower values than oxygen consumption time in the sediment (t c ). Figure 5 shows the dynamics during a wave cycle of the parameters u * , t D , DOU and

C w /C 0 for a 15 s period wave with U w =128 cm s -1 (u * m =5 cm s -1 , Re δ =2797
) and the highest oxygen consumption rate in the sediment (µ=2000 g m -3 d -1 , t c =432 s). The oxygen diffusion time through the DBL (t D ) varied from less than 2 s to about 40 s during each half a wave cycle, with t D =10 s.

SWI oxygenation is no longer limited by the transfer time through the DBL, resulting in a much higher DOU and interfacial concentration (DOU =2.1 g m -2 d -1 , C w =0.86 C 0 ), and a higher penetration depth (δ s =0.12 cm) compared to the laminar regime.

[ Finally, the penetration depth did not display any intra-wave dynamics. Cycle-averaged values for the oxygen concentration at the SWI also increase when the shear velocity increases (C w =0.78, 0.86, 0.91 C 0 for u * m =3, 5, 7 cm s -1 respectively). It is interesting to note that while the C w response to DOU fluctuations during a wave cycle is attenuated, its response to changes of cycle-averaged values of DOU is amplified.

Effect of the wave period and orbital velocity

[Figure 6] Fig. 7 shows the temporal evolutions of t D , DOU and C w /C 0 during a wave cycle for the same maximum shear velocity (u * m =5 cm s -1 ) and three wave periods (T =7, 12, 25 s), which correspond to three Reynolds numbers (Re δ =1627, 2404, 3865 resp.). For the three periods, the oxygen diffusion time (t D ) reached the same minimum value when the shear velocity is maximum, which occurred around the same phase in the three cases: the minimum DBL thickness is imposed by the maximum shear velocity which is the same for the three cases depicted on Fig. 7a. However, as the period increases, not only the time between two flow reversal increases, but also the laminar-turbulent transition occurs at an earlier phase during the accelerating phases of the wave cycle. Consequently, the duration of the turbulent periods (θ 1 + θ 2 ) is increased in absolute values and also in proportion of the wave period.

When the period increases, the duration of the relaminarization period θ 3 also increases from 2.2 to 5.1 s, which facilitates the relaminarization around flow reversal, yielding larger maximum values for t D . However, in proportion of the wave period, the duration of the period θ 3 decreases when the wave period increases. Regarding the DOU (Fig. 7b), when the wave period increases, the maximum value of the DOU increases while the minimum value decreases. The fluctuations of DOU and C w during the wave cycle (Fig. 7c) thus increase with the period. Although DOU cycle-averaged values vary little with the wave period (DOU ranging from 2.09 to 2.12 g m -2 d -1 ), DOU exhibits the minimum value for T =12 s. Cycle-averaged values for the oxygen concentration at the SWI display the same pattern with amplification (C w =0.86, 0.85, 0.88 C 0 for T =7, 12 and 25 s respectively).

[Figure 7]

Trends in cycle-averaged DO quantities

A similar dependence for DOU , VAR(DOU ), C w and VAR(C w ) on the wave characteristics was found for high (µ=2000 g m -3 d -1 ) and low (µ=50 g m -3 d -1 ) biogeochemical activities in the sediment.

The non-dimensional Sherwood number defined by Eq. ( 10) was used to describe the DO mass transfer at the SWI. It also provides a measure of the hydrodynamic efficiency of this mass transfer: the larger the Sherwood number, the lower the maximum shear velocity for the higher oxygen transfer at the SWI.

Since the Sherwood number was non-dependent on the biogeochemical activity of the sediment, it was used to summarize the dependency of the DO mass transfer at the SWI on the wave period and the orbital velocity (Fig. 8), and finally on the wave Reynolds number (Fig. ??). and 0.9 C 0 ). Over the range of wave conditions, Sh varied from 0.17 to 0.23, displaying minimum values for wave periods ranging from 10 to 15 s and orbital velocity ranging from 80 to 200 cm s -1 (Fig. 8a) corresponding to Reynolds numbers between 1500 and 3500 (Fig. ??a). The Sherwood number fluctuations during a wave cycle increased with the wave period and orbital velocity (i.e. with the Reynolds number), yielding VAR(Sh)=30 % (as a result of VAR(DOU )=17 % and VAR(C w )=1.75 %)

Laminar flow conditions (Re

for U w =200 cm s -1 and T =30 s (Re δ =6180).

[Figure 8] Fig. 9 shows DOU and C w /C 0 versus the cycle-averaged (mean) shear velocity u * during a wave cycle in the laminar and the turbulent regimes, and for two sediment biogeochemical activities. Whatever the flow regime, DOU is smaller and C w /C 0 is larger for the less active sediment (µ=50 g m -3 d -1 ) than for the more active sediment. In the laminar regime, DOU and C w /C 0 increased linearly with the mean shear velocity up to 1.6 cm s -1 . In the turbulent regime, both quantities reached an asymptote for the largest mean shear velocity.

[Figure 9]

As far as the oxygen penetration depth is concerned, in the laminar regime, it decreased as the wave period increased, whatever the orbital velocity and for the two biogeochemical activities (Fig. 10a). The oxic layer is thin, ranging from 0.38 to 0.79 cm and from 0.02 to 0.06 cm for µ=50 and µ=2000 g m -3 d -1

respectively. When the overlying flow is turbulent, oxygen penetration depth is thicker and again strongly depends on the biogeochemical activity in the sediment (δ s =1.03 and 0.15 cm for µ=50 and µ=2000 g m -3 d -1 respectively). However, in the turbulent regime, oxygen penetration depth no longer depends on the wave period and reaches an asymptote when shear velocity increases (Fig. 10b). For all the computations, DO penetration depth never showed intra-wave dynamics (i.e. VAR(δ s )=0).

[Figure 10] Finally, the convergence time (t st ) of DO profiles from an initial state (where C * =1 in water and C * =0 in sediment) to an equilibrium state was estimated from the calculations (data not shown). In the laminar regime, t st scaled from around 20 to 108 h and from around 20 to 92 h when biogeochemical consumption was low (µ=50 g m -3 d -1 ) and high (µ=2000 g m -3 d -1 ), respectively). Time to reach equilibrium state is shorter when the flow is turbulent (less than 15 h). For the less active sediment, t st varied little (between 11 and 13 h) with the wave period or orbital velocity. For the more active sediment, t st ranged from 1 to 15 h, and fastest convergence was obtained for the higher Reynolds number. These features were previously observed on a lake internal seiching with a period of 18 hours [START_REF] Lorke | Breathing sediments: The control of diffusive transport across the sediment-water interface by periodic boundary-layer turbulence[END_REF]. Below the SWI, fluctuations are rapidly damped out. In the laminar regime, however, no fluctuations of the DO concentration were observed, even in the overlying water.

Discussion

Three processes, all with different time scales, interplay during the wave period: (1) the oxygen consumption in the sediment, (2) the oxygen diffusion across the DBL, and (3) the DBL periodic dynamics. The relative importance of the first two processes was discussed by [START_REF] Higashino | Unsteady diffusional mass transfer at the sediment/water interface: Theory and significance for SOD measurement[END_REF] for steady flows. In a laminar regime, the oxygen diffusion time across the DBL (t D ) is much larger than the oxygen consumption time in the sediment (t c ): SWI oxygenation is limited by the oxygen supply, leading to low values for the DOU and C w . Conversely in a turbulent regime, the oxygen diffusion time across the DBL (t D ) is strongly reduced and SWI oxygenation is then scaled by the oxygen consumption in the sediment (t c ). In the case of unsteady wave forcing, the DBL periodicity whose time scale is the wave period (T ) will interfere with the oxygen consumption and diffusion.

For high-frequency waves with periods up to 30 s, the sediment oxygen consumption time (t c ) will always remain longer than the wave period (T ). Thus, any periodic fluctuation disappears rapidly in the sediment. In the turbulent regime, the oxygen diffusion time (t D ) can reach shorter values than the wave period (T ): every half a wave period, the t D value falls during the turbulent phases θ 1 and θ 2 .

Briefly put, the minimum value of t D should be compared to the duration of the turbulent part of the wave cycle (Fig. ??). which corresponds to sea swell conditions (10 < T < 15 s and 80 < U w < 200 cm s -1 ). In other words, some sea swells can be less efficient than wind waves (periods below 8 s) for oxygen mass transfer at the SWI.

As the wave

Comparison with DO dynamics under steady flows

Previous experimental studies investigated the DOU dynamics under steady flows and showed an increase of the uptake with the flow shear velocity [START_REF] Booij | Measuring the flux of oxygen to a muddy sediment with a cylindrical microcosm[END_REF][START_REF] Mackenthun | Effect of flow velocity on sediment oxygen demand: Experiments[END_REF][START_REF] Steinberger | Diffusional mass transfer at sediment-water interface[END_REF][START_REF] Tengberg | Intercalibration of benthic fluxes chambers I. Accuracy of flux measurements and influence of chamber hydrodynamics[END_REF]. These results were obtained for various sediment types (natural and artificial), but no indication of the sediment biogeochemical activities were given.

Consequently, direct comparison of our results with previous published experimental work cannot be quantitative. However, qualitatively, our study confirm an increase of the DOU with the shear velocity.

To our knowledge, the only study which can be used for a quantitative comparison is the numerical study by [START_REF] Higashino | Unsteady diffusional mass transfer at the sediment/water interface: Theory and significance for SOD measurement[END_REF]. DOU levels are quite similar, but the oxygenation at the interface is much greater in the case of steady flows. For u * =1 cm s -1 and µ=2000 g m -3 d -1 , oxygen transfer efficiency (i.e. Sherwood number) deduced from their calculations reaches Sh=0.36 whereas it reaches Sh=0.25 in our study.

High-frequency waves are less efficient than steady current for the oxygen mass transfer at the SWI. Such finding should meanwhile be debiased by the occurrence frequency of each flow condition.

Over the range of wave conditions we tested, the laminar regime corresponds to orbital velocities lower than 40 cm s -1 for the periods below 5 s and orbital velocities lower than 15 cm s -1 for periods up to 25 s. For the smallest periods corresponding to building wind waves (below 5 s), the wave height is also generally small and the laminar regime should be the most frequent regime in the bottom boundary layer, except in very shallow conditions (less than half a meter of water). This situation is likely to be encountered frequently in small closed basins (small lakes, lagoons, etc.).

As the wave period increases (when the wind fetch is larger), the laminar regime would become rarer: for instance, for a wave period of 8 s over 10 m water depth, the wave bottom boundary layer would be turbulent for a wave height H larger than 0.3 m, yielding Sh=0.18. For H=1.6 m, orbital velocities reach 60 cm s -1 (u * =2.1 cm s -1 , Sh=0.274), and yields on average a high SWI oxygenation (C w ≈ 0.8 C 0 ) due to large flux (DOU > 2 g m -2 d -1 ). Such wave conditions are frequently observed at sea (about 50 % of the time on the French Atlantic coast, [START_REF] Butel | Spatial variability of wave conditions on the French Atlantic coast using in-situ data[END_REF]. After [START_REF] Higashino | Unsteady diffusional mass transfer at the sediment/water interface: Theory and significance for SOD measurement[END_REF], reaching C w ≈ 0.8 C 0 under a steady flow requires a lower oxygen flux DOU =1.8 g m -2 d -1 associated with a lower shear velocity of 1.35 cm s -1 . Such steady shear velocity values correspond to a current speed of 30 cm s -1 at 1 m above a smooth bed. Even the most energetic tidal flow will not have such values more than 50 % of the time [START_REF] Butel | Spatial variability of wave conditions on the French Atlantic coast using in-situ data[END_REF]. Thus, it may be estimated that waves will promote a higher sediment oxygenation more frequently than steady currents. Additionally, DO profiles respond to the flow within a few hours. Even short events of wind waves or sea swell may notably stimulate shallow sediment oxygenation and oxygen exchange across the SWI.

Comparison with DO dynamics under low-frequency flows

Higashino et al. (2003) already attempted to describe the unsteady diffusional mass transfer under low-frequency periodic flows (periods ranging from 10 s to 10 h). Their numerical modelling was, however, restricted to the water column and they used a time varying water diffusivity derived from [START_REF] Dade | Near-bed turbulence and hydrodynamic control of diffusional mass transfer at the sea floor[END_REF] steady-state parameterization. This may be valid for slowly varying flows, but not for highfrequency waves. Time scales of diffusion across the viscous sublayer are of the same order of magnitude as the wave period for high-frequency flows. As a consequence the turbulent shear stresses display a phase lag across the boundary layer (Dohmen-Janssen, 1999) and the phase at which the laminarturbulent transition occurs during each acceleration of the wave cycle depends on the wave Reynolds number Re δ [START_REF] Jensen | Turbulent oscillatory boundary layers at high Reynolds numbers[END_REF]. These features are well reproduced by the turbulent boundary layer model we used to compute the unsteady diffusive forcing [START_REF] Guizien | 1DV bottom boundary layer modeling under combined wave and current: Turbulent separation and phase lag effects[END_REF].

To a small extent, the results presented in this paper tally with previous calculations of the oxygen distribution dynamics under low-frequency flows [START_REF] Higashino | Periodic diffusional mass transfer near sediment/water interface: Theory[END_REF], as intra-wave fluctuations of the Sherwood number effectively tends to zero when the period becomes shorter. Nonetheless, regarding the cycle-averaged value of the Sherwood number, our results are in contradiction to the afore-mentionned study. Authors reported a decrease of the cycle-averaged Sherwood number with the flow period, while we demonstrated a more complex pattern depending on the flow period and orbital velocity as well (Sh is minimum around T =10-15 s and U w =120-180 cm s -1 ). This difference is simply explained by the dependence of the maximum value of the Sherwood number on the wave period which was lacking in [START_REF] Higashino | Periodic diffusional mass transfer near sediment/water interface: Theory[END_REF] study. Only the intra-wave minimum value of the Sherwood number increased as the wave period decreased from 10 h to 10 s. Yet, we showed that the values of both the maximum and the minimum of the Sherwood number depend on the wave period (given a shear velocity maximum).

It should also be pointed out that our Sherwood number values (Sh around 0.2) are less than half of the values obtained by [START_REF] Higashino | Periodic diffusional mass transfer near sediment/water interface: Theory[END_REF] for similar flow conditions (Sh=0.5 for T =20 s,

u * m =1 cm s -1
). The authors extrapolated their results to state that when the wave period tends to zero, the Sherwood number should approach a maximum corresponding to the Sherwood number value for steady-state conditions at the shear velocity maximum. Using the simulations from a study of the DOU dynamics under steady flows [START_REF] Higashino | Unsteady diffusional mass transfer at the sediment/water interface: Theory and significance for SOD measurement[END_REF], the upper limit of the Sherwood number under periodic flow (when the maximum shear velocity is u * =1 cm s -1 ) should be Sh=0.36. We argue that Sh values were overestimated by [START_REF] Higashino | Periodic diffusional mass transfer near sediment/water interface: Theory[END_REF] probably because the interfacial oxygen concentration was assumed as a constant, while it was not in the present study.

Consequences for biological systems

Oxygen near the SWI is essential for both biological production and decomposition processes.

Oxygen concentration is a regulating mechanism for benthic community complexity [START_REF] Llansó | Effects of hypoxia on estuarine benthos: The lower Rappahannock River (Chesapeake Bay), a case study[END_REF], for regeneration of nutrients [START_REF] Rahm | On the mass transfer properties of the benthic boundary layer with an application to oxygen fluxes[END_REF], and it also constrains chemical [START_REF] Cai | Oxygen penetration depths and fluxes in marine sediments[END_REF] and biological [START_REF] House | Factors influencing the extent and development of the oxic zone in sediments[END_REF] reactions. In a laminar regime, the transfer of DO is strongly impeded by the thickness of the DBL and hypoxic conditions (C < 2 ml l -1 ≈ 2.4 g m -3 , Diaz and Rosenberg, 1995) may appear near the SWI. A continuous reduction of around 80-90 % in oxygen concentration at the SWI (C w =0.1-0.2 C 0 =1-2 g m -3 ) may significantly disturb the growth and metabolism of benthic organisms and can cause mass mortality of marine animals, resulting in benthic defaunation and fish community decline (Diaz and Rosenberg, 1995, and references herein). On the other hand, a turbulent regime promotes high levels of SWI oxygenation.

Oxygen concentration at the SWI can be dynamically driven by periodic fluctuations of the DBL thickness under wind waves and sea swell. Yet the amplitude of these fluctuations during a wave cycle may not be sufficient to stress the benthic macrofauna (around 1 % for C w /C 0 for instance). In the sediment, periodic fluctuations vanish rapidly and oxygen penetration depth largely depends on the oxygen consumption in the sediment and to a lesser extent, on the wave period in the laminar regime (Fig. 10a) and on the shear velocity in the turbulent regime (Fig. 10b). This sediment inertia is explained by the larger oxygen consumption time (t c ) in the sediment compared to the DBL dynamics time scale during a wave period (t w ). For flows having periods of hours (tides, inertial waves, etc.) or smaller consumption time (i.e. smaller C 0 , for instance), it can be expected that t w will reach values comparable to or even larger than t c . In such cases, periodic fluctuations of oxygen concentration may be expected not only in the overlying water but also in the sediment. However, the fluctuation range as well as the cycle-averaged value are not trivial to extrapolate since a strong coupling between the DBL dynamics, the oxygen diffusion across the DBL, and the oxygen consumption in the sediment is likely to exist when characteristic time scales (t D , t c and t w ) have similar orders of magnitude.

Finally, sediment oxygen consumption was modelled using a Michaelis-Menten kinetic formulation.

This formulation assumes explicitly that no acclimatization of the respiring community occurs during a wave cycle. Such absence of respiration acclimatization should be tested experimentally as a function of the wave period, for instance with high resolution DO measurements. Although promising results have accompanied the development of optical methods for oxygen concentration profiling [START_REF] Revsbech | A novel microsensor for determination of apparent diffusivity in sediments[END_REF], faster oxygen profiles measurements are still required to resolve the temporal dynamics of DO under wind waves and sea swell. Even the novel non-invasive eddy correlation technique [START_REF] Berg | Oxygen uptake by aquatic sediments measured with a novel non-invasive eddy-correlation techniqueroy[END_REF] requires an integration time of around 150 s, which is still too long to observe high-frequency fluctuations related to high-frequency wave turbulence. Turbulence and diffusion measurements should indeed be measured simultaneously with DO profiles at high-frequency to evaluate the tested model design and the biological interpretations of the mass transfer dynamics at the SWI under periodic waves.

Conclusions

The response of dissolved oxygen profiles to a wide range of wind waves to sea swell conditions over a smooth bottom was studied. Oxygen diffusion time across the diffusive boundary layer (t D ), dissolved oxygen uptake (DOU ), interfacial concentration (C w ), and oxygen penetration depth (δ s )

were described by their cycle-averaged value and their standard deviation during a wave cycle.

• As in steady flows, the flow regime (turbulent or laminar) controls both the DOU and the oxygen concentration at the SWI. In the laminar regime, oxygen mass transfer efficiency is strongly impeded, yielding small and constant DOU , C w and δ s values during the wave cycle. In the turbulent regime (frequent under waves), oxygen supply to the sediment is no longer limited resulting in higher values for DOU , C w /C 0 and δ s (Fig. 8, Fig. 10a andb).

• In the turbulent regime, DOU and C w exhibit periodic fluctuations in response to the unsteady DBL thickness: their dynamics were water controlled by the ratio between (a) the diffusion time across the DBL (t D ), and (b) the time allowed for this diffusion during a wave cycle (t w ). When t w (t) was smaller than t D (t), the DBL thickness varied more rapidly than the time required for diffusive flux to establish: diffusive flux corresponded to the average DBL thickness over a period. Conversely, when t w (t) was larger than t D (t), diffusive flux followed the local DBL thickness fluctuations (Fig. 5). Unsteady properties of the overlying high-frequency flow were never transmitted deeply into the sediment (VAR(δ s )=0).

• Dependence of the oxygen mass transfer efficiency on the wave period and the orbital velocity are summarized by the non-dimensional Sherwood number (Sh) trends. The latter was found nondependent on the sediment biogeochemical activity. The cycle-averaged Sherwood values (Sh) vary little, ranging from 0.17 to 0.23 and with minimum values for sea swells having periods ranging from 10 to 15 s (Fig. 8a). Intra-wave fluctuations VAR(Sh) increase as both wave period and orbital velocity (i.e. Reynolds number) increase, ranging from 0.1 % in the laminar regime to 30 % in the turbulent regime (Figs. 8b and ??b). Finally, although Sh values remained lower than values for a steady current with the same maximum shear velocity, high-frequency waves are likely to promote high sediment oxygenation more frequently than steady currents.

A Symbols

Symbol Unit Description C [M].[L] -3 Dissolved oxygen concentration C0 [M].[L] -3 Bulk water concentration C w [M].[L] -3 Interfacial concentration D [L] 2 .[T] -1 Vertical diffusivity DOU [M].[L] 2 .[T] -1 Dissolved oxygen uptake D s [L] 2 .[T] -1 Effective diffusivity Dt [L] 2 .[T] -1 Turbulent diffusivity D w [L] 2 .[T] -
[L] Oxygen penetration depth - Convergence criterion µ [M].[L] -3 .[T] -1 Biogeochemical activity of the sediment ν [L] 2 .[T] -1 Kinematic viscosity in water ν t [L] 2 .[T] -1 Eddy viscosity of the fluid φ - Sediment porosity ϕ [ o C] Wave phase ρ [M].[L] -3 Water density τ [M].[L] -1 .[T] -2 Bottom

  The new model was used to examine three fundamental questions about the dynamics of DO distribution under high-frequency waves. First, what are the characteristic time scales of the forcing that drive the intra-wave dynamics of DO ? Second, how does the DO mass transfer vary with the wave period and the orbital velocity (mean value and fluctuations during a wave cycle) ? And third, to what extent do wind waves and sea swell drive sediment oxygenation (DO penetration depth) ?

Fig. 4

 4 Fig. 4 displays computed turbulent diffusivity (i.e. eddy viscosity) profiles every 30 o along half a wave cycle in the turbulent regime. Temporal changes in D t values reflect the temporal variations during a wave cycle of the viscous sublayer thickness δ ν (t) (the area where D t < ν= 10 -2 cm 2 s -1 ).

  Figure 5]Values of t D exhibited periodic fluctuations every half a wave cycle (Fig.5b): it had a maximum around flow reversal and dropped rapidly at the onset of turbulence during flow acceleration. It decreased more gently until shear velocity reached a maximum (turbulent period labelled θ 1 ). As the shear velocity decreased, t D increased slightly during the turbulent period labelled θ 2 . Then it increased more rapidly as the bottom boundary layer became laminar during flow deceleration (relaminarization period labelled θ 3 ) to reach its maximum value again when the outer flow speed is zero. During θ 1 and θ 2 , t w is longer than t D (bold line on Fig.5b) while during θ 3 , t w is shorter than t D . While the t D minimum occurs when the shear velocity is at a maximum, the t D maximum occurs after the shear velocity minimum, close to the flow reversal point. This delay reflects the phase lagging of turbulent diffusivity across the bottom boundary layer. As a consequence of the t D fluctuations, DO concentrations at the SWI also displayed intra-wave dynamics: the DOU at the SWI varied between a minimum and maximum value during each half wave cycle, yielding VAR(DOU )=14 % (Fig.5c). The minimum value of the DOU occurs at the end of the relaminarization period θ 3 . The DOU starts to increase at the laminar-turbulent transition when diffusion time had strongly decreased, reaches a maximum during the turbulent period, and then decreased until the next laminar-turbulent transition. It is noteworthy that during θ 3 , the DOU decreases steadily without reflecting the large fluctuations of t D around flow reversal. As a result of the DOU fluctuations, oxygen concentration at the SWI exhibited similar intra-wave fluctuations, although those were damped (VAR(C w /C 0 )=1 %, Fig.5d) and delayed compared to the DOU fluctuations.

Fig. 6

 6 Fig.6shows the temporal evolutions of t D , DOU and C w /C 0 during a wave cycle for the same wave period (T =15 s) and three orbital velocities (U w =62, 128, 195 cm s -1 ), which correspond to three maximum shear velocities (u * m =3, 5, 7 cm s -1 respectively) and three Reynolds numbers (Re δ =1355, 2797, 4261 resp.). When the maximum shear velocity increases, the minimum value for t D decreases from 5.3 to 0.7 s since the minimum DBL thickness is imposed by the maximum shear velocity. Moreover, the onset of turbulence occurs at an earlier phase during the accelerating phases of the wave cycle and the duration of period θ 1 increases. The duration of the relaminarization period θ 3 decreases from 5 to 2.6 s which limits the laminar DBL development around flow reversal. The maximum values for t D around flow reversal decreases, although remaining larger than the wave period for the three cases

  δ ≤ 320 on Fig. ?? and light grey areas on Fig. 8) significantly hinder the transfer of dissolved oxygen and strongly reduce the presence of oxygen at the SWI, yielding very low non-dimensional mass transfer coefficients on average (Sh < 0.05). Intra-wave cycle fluctuations are negligible (VAR(Sh) < 0.1 %, Figs. 8b and ??b). When the wave boundary layer becomes turbulent (Re δ > 320), DOU and C w displayed larger values than in the laminar regime (up to 2.15 g m -2 d -1

4. 1

 1 Characteristic time scales of non-steady DO dynamics under oscillatory flows Our simulations showed that in the turbulent regime, the periodic fluctuations of the DBL thickness under oscillatory forcing induced periodic fluctuations of the DO concentration above and at the SWI. These fluctuations are not phase locked, evidence of transient adjustment to DBL thickness fluctuations.

  Reynolds number increases, the minimum value of t D decreases and the duration of the turbulent phases increases (when t w > t D ). DOU has then more time to adjust to the lower values of the minimum of t D , and its intra-wave fluctuations grow. However, when the minimum value of t D remains longer than the duration of the turbulent part of the wave cycle ((θ 1 + θ 2 )/min(t D ), 2050 < Re δ < 2500 on Fig. ??), the maximum value of the DOU will not have time to equilibrate with the low value of t D . The maximum values of the Sherwood number are limited by the short duration of the turbulent phases. Consequently, intra-wave fluctuations are supported by the minimum values of the Sherwood number which decrease rapidly. Cycle-averaged value of the Sherwood number (Sh) decreases. As the wave Reynolds number further increases, the minimum value of t D decreases until reaching shorter values than the period θ 1 (θ 1 /min(t D ), 2800 < Re δ < 3300 on Fig. ??). The maximum value of the Sherwood number increases rapidly, following the decrease of the t D minimum, while the minimum value of the Sherwood number is limited by the duration of the relaminarization period θ 3 . Cycle-averaged value of the Sherwood number (Sh) increases. As a result of the competition between the oxygen diffusion time and the time granted to this rapid diffusion, Sh exhibits a minimum value for Re δ between 2050 and 3300 (light gray area on Fig. ??a),

Figure 1 :)Figure 2 :

 12 Figure 1: Schematic representation of non-dimensional dissolved oxygen concentration profile (C * =C/C0) at the sediment-water interface (SWI).

Figure 3 :

 3 Figure 3: Maximum shear velocity u * m dependence on orbital velocity Uw and period T of the wave. Dark grey delimits the unstable breaking waves (H < 0.8 D) and light grey delimits the laminar regime (Re δ < 320).

Figure 4 :Figure 5 :Figure 6 :Figure 7 :

 4567 Figure 4: Vertical profiles of turbulent diffusivity D t at different wave phases ϕ, for u * m =5 cm s -1 and T =25 s.

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Isograms of (a) Sh, and (b) VAR(Sh) showing dependences on orbital velocity Uw and wave period T . Dark grey delimits the unstable breaking waves (H < 0.8 D) and light grey delimits the laminar regime (Re δ < 320).

Aknowledgements

We would like to thank Jennifer Guarini for correcting the English and for commenting on an earlier version of this paper. The critical remarks of two anonymous reviewers helped to improve the manuscript. This work is a partial fulfillment of the PhD Thesis of Mathieu Chatelain who was supported by a doctoral fellowship from the French National Center for Scientific Research (CNRS).

List of Figures


domain: sdv sde

text:
Catharina J M Philippart
email: katja.philippart@nioz.nl

Ana Amaral

Ragnhild Asmus

Judith Van Bleijswijk

Julie Bremner

Fred Buchholz

Miguel Cabanellas-Reboredo

Diana Catarino

Andre Cattrijsse

François Charles

André Cattrijsse

Thierry Comtet

Alexandra Cunha

Salud Deudero

Jean-Claude Duchêne

Simonetta Fraschetti

Franck Gentil

Arjan Gittenberger

Katell Guizien

João M Gonçalves

Giuseppe Guarnieri

Iris Hendriks

Birgit Hussel

Raquel Pinheiro Vieira

Bastian T Reijnen

Iris Sampaio

Ester Serrao

Isabel Sousa Pinto

Eric Thiebaut

Frédérique Viard

Alain F Zuur

Spatial synchronies in the seasonal occurrence of larvae of oysters (Crassostrea gigas) and mussels (Mytilus edulis/galloprovincialis) in European coastal waters

Keywords: spatial synchrony, synoptic scale, larval occurrence, Crassostrea gigas, Mytilus edulis, Mytilus galloprovincialis

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

Evolution in marine invertebrates with complex life-cycles, such as bivalves, has resulted in synchronised reproductive cycles within species. The ultimate selective advantage of synchronization is to produce larvae under those circumstances which maximise their survival [START_REF] Olive | Annual breeding cycles in marine invertebrates and environmental temperature: Probing the proximate and ultimate causes of reproductive synchrony[END_REF]. Knowledge of temporal patterns of occurrence of larval stages of various species under different ecological settings may help, therefore, to identify environmental factors that enhance species-specific larval survival. This information can only be gained from robust comparison of heterogeneity across broad geographical scales, and is crucial for predicting the consequences for invertebrate populations under environmental changes such as eutrophication and global climate change [START_REF] Olive | Annual breeding cycles in marine invertebrates and environmental temperature: Probing the proximate and ultimate causes of reproductive synchrony[END_REF][START_REF] Harley | The impacts of climate change in coastal marine systems[END_REF][START_REF] Goberville | Climate-driven changes in coastal marine systems of western Europe[END_REF]. Synchronized reproduction has often been linked with water temperature, which varies with season and latitude [START_REF] Orton | Sea-temperature, breeding and distribution in marine animals[END_REF][START_REF] Giese | Comparative physiology -annual reproductive cycles of marine invertebrates[END_REF]. However, seasonal cycles of photoperiod or changes in food availability may also steer the timing of reproduction of marine benthic invertebrates [START_REF] Starr | Direct coupling of marine invertebrate spawning with phytoplankton blooms[END_REF][START_REF] Bentley | Physiological and environmental control of reproduction in Polychaetes[END_REF]. Species with latitudinally separated populations may have, therefore, predictably different reproductive cycles [START_REF] Orton | Sea-temperature, breeding and distribution in marine animals[END_REF][START_REF] Thorson | Reproduction and larval development of Danish marine bottom invertebrates[END_REF][START_REF] Dekshenieks | Environmental effects on the growth and development of eastern oyster, Crassostrea virginica (Gmelin, 1791), larvae -a modeling study[END_REF][START_REF] Hofmann | Modeling oyster populations II. Adult size and reproductive effort[END_REF][START_REF] Bhaud | Some examples of the contribution of planktonic larval stages to the biology and ecology of polychaetes[END_REF].

If waters are warm and food availability is relatively constant, e.g. in the tropics, bivalves generally display multiple spawning events within one year [START_REF] Pouvreau | Application of a dynamic energy budget model to the Pacific oyster, Crassostrea gigas, reared under various environmental conditions[END_REF]. In continuously cold waters where food availability is strongly seasonal, such as polar regions, bivalves and other invertebrates often show one strong spawning peak per year (Thorson, 1959;[START_REF] Kang | The annual reproductive pattern of the Antarctic clam, Laternula elliptica from Marian Cove, King George Island[END_REF]. In temperate waters where both temperature and food availability are seasonal, the timing of spawning of bivalves is also strongly seasonal.

Furthermore, boreal-temperate invertebrates tend to spawn earlier and for a more extended period in the warmer part than in the colder part of their distribution range (Thorson, 1959;[START_REF] Seed | Ecology[END_REF][START_REF] Bauer | Testing generalizations about latitudinal variation in reproduction and recruitment patterns with sicyoniid and caridean shrimp species[END_REF][START_REF] Defeo | Macroecology of population dynamics and life history traits of the mole crab Emerita brasiliensis in Atlantic sandy beaches of South America[END_REF][START_REF] O'riordan | Spatial variation in the recruitment of the intertidal barnacles Chthamalus montagui Southward and Chthamalus stellatus (Poli) (Crustacea: Cirripedia) over an European scale[END_REF].

Many of the findings corroborating a latitudinal gradient in environmental conditions and in bivalve reproduction come from large-scale (100s-1000s of km) observations in marine waters along long and uninterrupted north-south orientated coastlines such as the Pacific coast of North America [START_REF] Bauer | Testing generalizations about latitudinal variation in reproduction and recruitment patterns with sicyoniid and caridean shrimp species[END_REF] and the Atlantic coast of South America [START_REF] Defeo | Macroecology of population dynamics and life history traits of the mole crab Emerita brasiliensis in Atlantic sandy beaches of South America[END_REF].

More detailed studies along these coastlines, however, often revealed mesoscale (10s-100s of km) variation in reproductive factors which appears to be related to coastal topographic features such as bays and headlands (e.g., [START_REF] Lester | Reproduction on the edge: large-scale patterns of individual performance in a marine invertebrate[END_REF]. These features structure local hydrodynamics such as residence times [START_REF] Frantzen | Recruitment of blue mussels, Mytilus edulis L., on suspended collectors in Finnmark, North Norway (70-71°N)[END_REF] and upwelling events (e.g., [START_REF] Brink | The near-surface dynamics of coastal upwelling[END_REF][START_REF] Ebert | Latitudinal variation in size structure of the west coast purple sea urchin: a correlation between headlands[END_REF][START_REF] Connolly | A latitudinal gradient in recruitment of intertidal invertebrates in the northeast Pacific ocean[END_REF][START_REF] Largier | Considerations in estimating larval dispersal distances from oceanographic data[END_REF]. The subsequent local and mesoscale variation in the timing of bivalve reproduction and in larval distribution may interfere with synoptic-scale relationships between latitude, environmental conditions and the seasonality in the presence of larvae in coastal waters.

In contrast to the Pacific coast of North America and the Atlantic coast of South America, Europe does not have an uninterrupted north-south orientated coastline because it comprises semi (e.g. North Sea) to almost fully (e.g. Baltic Sea, Mediterranean and Black Sea) land-locked seas. This topographical complexity may result in a decoupling of latitude and environmental conditions. Chlorophyll-a concentrations tend to be higher in coastal than open seas as the result of nutrient supply by rivers or upwelling events (e.g., [START_REF] Babin | Variations in the light absorption coefficients of phytoplankton, nonalgal particles, and dissolved organic matter in coastal waters around Europe[END_REF]. The enclosed seas have different SSTs compared to the more open marine waters of Europe at similar latitudes [START_REF] Belkin | Rapid warming of large marine ecosystems[END_REF].

In this paper, we investigated the presence of bivalve larvae at 12 sites in European coastal waters, ranging from 37°N to 60°N (distance ≈ 2400 km) and from 27°W to 18°E (distance ≈ 4000 km), during one year of observations. Although previous studies have examined mesoscale (100s of km) relationships between latitude, environmental conditions and invertebrate reproduction in European marine waters (e.g., [START_REF] Haury | Patterns and processes in the time-space scale of plankton distribution[END_REF][START_REF] Brown | Geographical variations in the reproduction of the horse mussel, Modiolus modiolus (Mollusca: Bivalvia)[END_REF][START_REF] Lewis | Latitudinal trends in reproduction, recruitment and population characteristics of some rocky littoral molluscs and cirripedes[END_REF][START_REF] Jenkins | Spatial and temporal variation in settlement and recruitment of the intertidal barnacle Semibalanus balanoides (L.) (Crustacea: Cirripedia) over a European scale[END_REF][START_REF] O'riordan | Spatial variation in the recruitment of the intertidal barnacles Chthamalus montagui Southward and Chthamalus stellatus (Poli) (Crustacea: Cirripedia) over an European scale[END_REF][START_REF] Cardoso | Spatial variability in growth and reproduction of the Pacific oyster Crassostrea gigas (Thurnberg, 1793) along the west European coast[END_REF][START_REF] Burrows | Spatial variation in size and density of adult and post-settlement Semibalanus balanoides: effects of oceanographic and local conditions[END_REF], our study was the first which examined the seasonality in bivalve larvae at a synoptic scale (1000s of km) in this part of the world.

We assumed that the presence of larvae in the water is a proxy for the timing of reproduction. We analysed the larval samples for the presence of two bivalve species, i.e. Pacific oysters (Crassostrea gigas) and mussels (Mytilus edulis/galloprovincialis). We tested how much of the observed seasonal patterns of presence of bivalve larvae at 12 European stations could be explained by various models, including one that stated that there was no synchronisation in seasonality of larval presence at all between the locations (null hypothesis), one that assumed that there was one common seasonality pattern for all stations within Europe, and various models which supposed that the variation in seasonality could be grouped according to specific spatial scales, taxonomic groups, or several combinations of these factors.

Materials and methods

Sampling

Sampling was carried out at 12 coastal sites in Europe. Sites were selected to maximise the range of environmental conditions and latitude in European coastal waters (Table 1). All study sites were permanently submerged, and preferably far from freshwater inputs and other major sources of impacts such as dredging activities and harbours. In addition, the sites were accessible throughout the year to be able to sample at the required frequency and for the full study period. Sites consisted of a single station except for the Azores where two stations were sampled.

To standardize sampling and laboratory procedures, the Royal Netherlands Institute for Sea Research (NIOZ) supplied all materials to filter and store the samples. Sampling started on December 15, 2006, and was aimed to take place once every two weeks throughout the year until December, 2007. Water samples of 10 L were taken from the water surface using a bucket, preferably at high tide for tidal stations and around noon for non-tidal stations. One subsample of 2 L was taken from the bucket, while the water was gently stirred to keep the larvae into suspension. This subsample was filtered over a 31.5 µm-mesh nylon filter with a diameter of 47 mm (tailormade by "Merrem & la Porte BV", Zaltbommel, The Netherlands) using minor suction of 0.5 to 1 bar. Material was washed down onto the filter with 0.2 µmfiltered seawater with salinity comparable to that of the samples. The filters were stored separately in pre-labelled cryovials at -80°C. At the end of the field campaign, the frozen samples were transported by couriers on dry-ice to NIOZ and immediately stored at -80 ⁰C after arrival.

For several stations, we measured chlorophyll-a concentrations and sea surface temperature (SST) at or close to the larvae sampling stations (i.e., Sylt, Texel, Oostende and Lecce). Data on environmental conditions for Roscoff and Banyuls were provided by the French Coastal Environmental Monitoring Service (SOMLIT). For all other stations, monthly water temperatures (SST) were derived from monthly averaged satellite data (http://oceancolor.gsfc.nasa.gov) from December 2006 to December 2007.

Species

In the mid-1960s, the Pacific oyster (C. gigas) was introduced for commercial purposes from the Northwest Pacific to several locations in Europe. Since then it has spread widely and wild populations can now be found along most of Europe's coastlines, from the northern Mediterranean to southern Norway [START_REF] Miossec | Crassostrea gigas (Pacific oyster)[END_REF]. In northern European waters, there has been a transition from triploid sterile oyster populations sustained by hatchery production to recruitment from wild-settlement [START_REF] Gollasch | Introduced aquatic species of the North Sea coasts and adjacent brackish waters[END_REF]).

The absence of Pacific oysters in northern European waters may be caused by too low temperatures for survival and reproduction. Pacific oysters cannot survive at temperatures below 3°C [START_REF] Kobayashi | A population dynamics model for the Japanese oyster, Crassostrea gigas[END_REF][START_REF] Cardoso | Spatial variability in growth and reproduction of the Pacific oyster Crassostrea gigas (Thurnberg, 1793) along the west European coast[END_REF]. While Pacific oysters initiate gametogenesis between 8 to 11 °C, they spawn after a threshold temperature of 18 to 20 °C has been reached [START_REF] Mann | Some biochemical and physiological aspects of growth and gametogenesis in Crassostrea gigas and Ostrea edulis grown at sustained elevated temperatures[END_REF][START_REF] Fabioux | Temperature and photoperiod drive Crassostrea gigas reproductive internal clock[END_REF]. For optimal development, oyster larvae need water temperatures of at least 22°C [START_REF] Kobayashi | A population dynamics model for the Japanese oyster, Crassostrea gigas[END_REF][START_REF] Shatkin | Considerations regarding the possible introduction of the Pacific oyster (Crassostrea gigas) to the Gulf of Maine: A review of global experience[END_REF][START_REF] Rico-Villa | A flow-through rearing system for ecophysiological studies of Pacific oyster Crassostrea gigas larvae[END_REF]. Reproduction appears to be more sensitive to low temperatures than survival. For example, although Pacific oysters can grow such as the Bay of Morlaix (Roscoff sampling station), they are not able to reproduce in these relatively cold waters [START_REF] Lejart | Etude du processus invasif de Crassostrea gigas en Bretagne: état des lieux, dynamique et conséquences écologiques[END_REF].

The southern edge of distribution may also be temperature-limited. The temperature tolerance of Pacific oysters ranges between 3°C and 35°C [START_REF] Cardoso | Spatial variability in growth and reproduction of the Pacific oyster Crassostrea gigas (Thurnberg, 1793) along the west European coast[END_REF], and summer adult mass mortalities have been documented at hatcheries in France and Ireland since the mid-1980s and mid-1990s, respectively [START_REF] Cotter | Summer mortality of the Pacific oyster, Crassostrea gigas, in the Irish Sea: The influence of growth, biochemistry and gametogenesis[END_REF][START_REF] Huvet | Reproductive effort of Pacific oysters: A trait associated with susceptibility to summer mortality[END_REF].

The Mytilus edulis complex comprises three nominal species, namely M. edulis, M. galloprovincialis and M. trossulus (Śmietanka et al., 2009). Whilst M. trossulus inhabits the Baltic Sea and Loch Etive in Scotland, M. edulis occurs from the Barents Sea in the north to the Atlantic coast in France to the south, and M. galloprovincialis extends from the Orkney Island in the north to the Mediterranean and Black Sea in the south [START_REF] Gosling | The mussel Mytilus: ecology, physiology genetics and culture[END_REF][START_REF] Berge | Ocean temperature oscillations enable reappearance of blue mussels Mytilus edulis in Svalbard after a 1000 year absence[END_REF][START_REF] Zbawicka | Scottish Mytilus trossulus mussels retain ancestral mitochondrial DNA: complete sequences of male and female mtDNA genomes[END_REF]. These mussel species can also interbreed and produce fertile hybrids where their distributions overlap, i.e. the French Atlantic coast up to northern Scotland [START_REF] Bierne | Introgression patterns in the mosaic hybrid zone between Mytilus edulis and M. galloprovincialis[END_REF][START_REF] Śmietanka | Mitochondrial DNA lineages in the European populations of mussels (Mytilus spp.)[END_REF], 2009). Whilst M. edulis is genetically homogeneous throughout its range, M.

galloprovincialis is subdivided into a Mediterranean group and an Atlantic group with a break point at the Almeria-Oran oceanographic front in the Mediterranean Sea [START_REF] Quesada | A multilocus allozyme discontinuity in the mussel Mytilus galloprovincialis: the interaction of ecological and life-history factors[END_REF]. M. galloprovincialis is slowly spreading northwards invading territory once exclusive to M. edulis which may be a result of global climate change [START_REF] Beaumont | Blue Mussel Mytilus edulis, Mediterranean mussel M. galloprovincialis[END_REF]. Adults of Mytilus species are virtually absent in the Azores [START_REF] Hawkins | Thoughts on the ecology and evolution of the intertidal biota of the Azores and other Atlantic islands[END_REF] and whether larvae of mussels are present is still unknown.

Molecular analyses

Total DNA was extracted by the UltraClean™ Soil DNA Isolation Kit (MOBIO Laboratories Inc.) using the alternative protocol for maximum yields (version 03252005). The filter was incised by a sterile scissor to increase the surface area before being subjected to the steps in the extraction protocol.

The DNA extracts were checked on a 1% agarose gel stained with ethidium bromide and visualized. The primers did not match with the male lineage of M. edulis/galloprovincialis (two mismatches in the final four bases at the 3' side). For M. trossulus (both females and males), there were no mismatches in the final four bases at the 3' side, but three mismatches further on.

Specific primers for

Each PCR reaction mixture contained: 2.5 μL BioTherm™ PCR Buffer (contains 15mM MgCl2), 2.5 μL DNTP's (Bioline, 2.5 mM), 0.25 μL of each primer (50 μM), 0.5 μL BSA (20 mg ml -1 ), 0.125 μL BioThermPlus™ Taq (5 units μl -1 ), 15.875 μL PCR H 2 O, 3 μL DNA template. After an initial denaturation step at 94°C for 120 seconds, the PCR programmes followed 36 cycles at 94°C for 30 seconds, annealing at 52.5°C (C. gigas) or at 50°C (M. edulis/galloprovincialis) for 30 seconds, extension at 72°C for 60 seconds and a final extension step at 72°C for 420 sec.

Presence of PCR inhibitors in a sample (e.g. from the seawater) may result in a false negative, i.e. a negative amplification result whilst larvae were present in that sample. Our results were therefore checked by means of PCR inhibition tests. For this test we used the protocol as described above, but instead of 3 μl sample we used 2 μl of the diluted sample (1:1, 1:5, 1:25 and 1:125) and added 1 μl of positive control consisting of adult DNA of each target species and visualized the results on gel.

The presence/absence of larvae of the target species was scored from the gels immediately after PCR amplification. Presence was detected if there was at least one D-larva (early veliger larval stage of bivalves) larger than 50 μm in the 2-L sample. This implied that it was possible to score larvae from mean concentrations of 0.5 -2.5 larvae per litre (= 500 -2.500 larvae per m 3 ) onwards. We scored the visibility of the bands as an index of presence: larvae were assumed to be absent (0) if no band was visible, and present (1) if the bands were weakly to strongly visible. Molecular analyses, including scoring, were performed by only one person (BTR). and the Mediterranean (which is subdivided into 7 ecoregions).

Spatial scales

Following the protocol in [START_REF] Zuur | A protocol for data exploration to avoiding common statistical problems[END_REF], an initial data exploration was applied. Hereafter, we modelled the data by means of a generalised additive mixed effects model (GAMM) with a binomial distribution, i.e.:

( ,1) logit( ) covariates is is is i is M Binomial p p a       
where M is is the observed presence or absence of Crassostrea gigas or Mytilus edulis/galloprovincialis at site i at time s, and p is is the probability that larvae of oysters or mussels are present at site i at time s. The random intercept a i allows for a random variation around the intercept α, and is assumed to be normally distributed with mean 0 and variance σ 2 site . The random intercept ensures that observations from the same site are allowed to be correlated. Due to the short length of the time series, we did not apply more advanced correlation structures like auto-regressive correlation [START_REF] Zuur | Mixed effects models and extensions in ecology with R[END_REF], but instead only used the symmetric compound correlation as imposed by the random intercept.

Based on the underlying questions (see section 2.4), a series of hypotheses and models were used for the 'covariates' part (Table 3). The smoothers f(.) were estimated using the R package (R Development Core Team, 2009), gamm4 [START_REF] Wood | gamm4: Generalized additive mixed models using mgcv and lme4. R package version 0[END_REF] and cubic regression splines. The covariate Subspecies was fitted as a categorical variable. In order to fit all models on the same data, the observations from the Azores (Faial North and South) were removed from the analysis.

To compare models, we followed an Information Theoretic approach [START_REF] Anderson | Avoiding pitfalls when using information-theoretic methods[END_REF]. We calculated differences Δ i between the Akaike Information Criterium (AIC) of each model and the minimum AIC. [START_REF] Anderson | Avoiding pitfalls when using information-theoretic methods[END_REF] state that the level of empirical support for model i is:

 substantial if Δ i is between 0 and 2 (hence, these are models with similar AICs as the optimal model),  considerable less if Δ i is between 4 and 7, and

 essentially none if Δ i is larger than 10.
Akaike weights w i were calculated following [START_REF] Anderson | Avoiding pitfalls when using information-theoretic methods[END_REF], which place the Δ i on a more interpretable scale by using:

1 2 1 2 i i i i e w e      
These weights have the convenient ability that they can be interpreted as probabilities that a given model is judged the best model on repeated sampling. If the weight for a particular model has a value of 0.75, for example, this implies that this model has a probability of 75% of being the best model within the series tested. Note that the AIC values obtained by a GAMM with a binomial distribution should be considered as approximate.

Results

Oysters (Crassostrea gigas)

Larvae of C. gigas were observed at 5 of the 12 sites, and occurred more often at Sylt, Texel and Roscoff than at Oostende and Algarve (Fig. 2).

The values of the Akaike weights w i indicated that model H1 has the highest probability (58%) of being the best model within the series tested, followed by model H2 (25%) and model H3 (17%) (Table 3).

The best model explaining the presence/absence of oyster larvae in European coastal waters was the one that assumes one common seasonal trend at all stations (H1 in Table 4). This model suggested that oyster larvae are most likely to occur in summer with peak values around early August (day number experiment ≈ 230) (Table 5; Fig. 3).

The other model for which the level of empirical support is substantial (i.e. the value of Δi is between 0 and 2, so its AIC is considered to be similar to that of the optimal model H1) was the one that assumes one underlying seasonal trend and one latitudinal effect for all the time series (H2; Table 5; Fig. 4) with higher probability of occurrence in the north than in the south.

The model that assumes that the latitude effect changes over time (H3) had some empirical support as well (Δi = 2.4). The results of model H3 suggested that larvae have a higher probability of occurrence at higher latitudes and exhibit a seasonal pattern for day number with earlier occurrences at higher latitudes in the range of those tested (Fig. 5).

Mussels (Mytilus edulis/galloprovincialis)

Mussel larvae were observed at all 12 sites, and were present during most of the year at Povoa, Sylt, Texel, Oslo and Banyuls but were rarely observed at Helgoland and Mallorca (Fig. 2). The values of the Akaike weights w i indicate that that model H4 has the highest probability (52%) of being the best model within the series tested, followed by model H5 (25%)

(Table 3).

The best model explaining the presence/absence of mussel larvae in European coastal waters was the one that assumed four different seasonal trends, one for each LME (H4 in Table 3). For the North Sea and the Iberian Coast, this model suggested that the highest probability of occurrence of mussel larvae was in late winter (end of February) and around the end of September (Table 6; Fig. 6). The Celtic-Biscay Shelf was characterised by an increase in probability of occurrence of mussel larvae between February and June, whilst the Mediterranean shows a decrease in probability of occurrence of mussel larvae throughout spring and summer (Fig. 6).

The other model for which the level of empirical support is substantial (i.e. the value of Δi is between 0 and 2) is the one that assumes that all the time series have one underlying seasonal trend (H1 in Table 3) with peak values around mid-April (day number experiment ≈ 120) (Table 6; Fig. 7).

Discussion

At optimal water temperatures and feeding conditions, larvae of oysters and mussels may reach metamorphosis within two weeks [START_REF] Satuito | Larval development of the mussel Mytilus-edulis galloprovincialis cultured under laboratory conditions[END_REF][START_REF] Rico-Villa | Influence of food density and temperature on ingestion, growth and settlement of Pacific oyster larvae, Crassostrea gigas[END_REF]. This implies that only under such maximum growth conditions for larvae, reproductive periods could have been missed for those stations were we were able to sample as planned, i.e. once every fortnight.

Due to the presence of a detection limit of larvae in our samples, we cannot distinguish between the presence of larvae at low densities (< 500 larvae per m 3 ) and the actual absence across sampling sites and/or sampling periods. We might, therefore, have missed the presence of larvae if their densities were low, e.g. as the result of low reproductive output of bivalve stocks and/or as the result of dilution if the sampling site was not close to the nearest spawning stock.

In addition, our sampling and analyses techniques were not quantitative and did not allow us to discriminate between a few larvae exceeding the threshold or higher larval densities. Interpretation of the results should, therefore, be done within the light of these limitations.

The strength of our methods, however, lies within the standardisation of the sampling techniques throughout Europe and the consistency of the analyses of the samples which were performed in one laboratory by one person. This reduced the possible noise in the data as the result of variation in sampling and analysis techniques between laboratories to a large extent.

Oysters (Crassostrea gigas)

The best model on the occurrence of oyster larvae in European coastal waters (Fig. 3) suggested that the reproduction of oysters most probably exhibit one common seasonal trend at all stations where larvae were observed, with highest probability of larval occurrence around early August.

This is in agreement with previous observations that natural populations of Pacific oysters in north-western Europe generally spawn in summer and autumn [START_REF] Pouvreau | Application of a dynamic energy budget model to the Pacific oyster, Crassostrea gigas, reared under various environmental conditions[END_REF][START_REF] Cardoso | Spatial variability in growth and reproduction of the Pacific oyster Crassostrea gigas (Thurnberg, 1793) along the west European coast[END_REF][START_REF] Enríquez-Díaz | Gametogenesis, reproductive investment, and spawning behavior of the Pacific giant oyster Crassostrea gigas: evidence of an environmentdependent strategy[END_REF][START_REF] Troost | Effects of an increasing filter feeder stock on larval abundance in the Oosterschelde estuary (SW Netherlands)[END_REF][START_REF] Dutertre | Rising water temperatures, reproduction and recruitment of an invasive oyster, Crassostrea gigas, on the French Atlantic coast[END_REF].

The second best but comparably good model (Fig. 4) suggested that larvae exhibit one seasonal trend and one underlying latitudinal effect in probability of occurrence, while the third best model (Fig. 5) suggested that oyster larvae also occurred earlier at higher latitudes in the range of those tested. If the timing of reproduction is predominantly triggered by a threshold in temperature, it would be expected that spawning would occur earlier in the south than in the north. Our findings were in agreement with previous observations by Cardoso and co-authors (2007), however, who showed that spawning events (identified by sudden decreases in gonadal mass) of Pacific oysters occurred earlier in northern than in southern sites along the west European coast (from 46 to 53°N).

The above evidence suggests that the onset of spawning is not only determined by temperature but also by other factors, possibly food availability.

Energy flow models of Crassostrea virginica indicated, for example, that both temperature and food supply affect the timing of the spawning season and the reproductive output [START_REF] Dekshenieks | Environmental effects on the growth and development of eastern oyster, Crassostrea virginica (Gmelin, 1791), larvae -a modeling study[END_REF][START_REF] Hofmann | Modeling oyster populations II. Adult size and reproductive effort[END_REF]. Within

NW-Europe, increased food conditions in the north were suggested to advance the timing of spawning of Pacific oysters during the warmest months of the year [START_REF] Cardoso | Spatial variability in growth and reproduction of the Pacific oyster Crassostrea gigas (Thurnberg, 1793) along the west European coast[END_REF]. Their suggestion on the importance of food availability on the timing of spawning of Pacific oysters were corroborated by our findings that chlorophyll a concentrations were considerably higher at the northern locations than at Roscoff (Fig. 2).

Larvae of Crassostrea gigas were detected by molecular techniques in winter and spring at the sampling site near Sylt (Fig. 2) which was confirmed by microscopic analyses of water samples taken at the same time (R. Asmus, unpublished data). These observations suggest that spawning can occur at much lower temperatures (less than 10⁰C; Fig. 2) than previously indicated as threshold values [START_REF] Mann | Some biochemical and physiological aspects of growth and gametogenesis in Crassostrea gigas and Ostrea edulis grown at sustained elevated temperatures[END_REF][START_REF] Fabioux | Temperature and photoperiod drive Crassostrea gigas reproductive internal clock[END_REF][START_REF] Lejart | Etude du processus invasif de Crassostrea gigas en Bretagne: état des lieux, dynamique et conséquences écologiques[END_REF]. Whether this might be the result of rapid adaption of the invasive Pacific oysters towards a relatively new environment [START_REF] Prentis | Adaptive evolution in invasive species[END_REF] or that Sylt harbours a different subpopulation of Pacific oysters remains to be investigated, e.g. by determining the traits with regard to phenology and larval survival of oysters from various origins and possible subpopulations.

Mussels (Mytilus edulis/galloprovincialis)

The timing of reproduction and the duration of the reproductive cycle in

Mytilus edulis, M. galloprovincialis and their hybrids have been shown to vary both spatially and temporally under local environmental conditions like temperature and food availability [START_REF] Wilson | Reproduction in Mytilus edulis L. (Mollusca: Bivalvia) in Carlingford Lough, Northern Ireland[END_REF][START_REF] Bayne | Reproduction in bivalve molluscs under environmental stress[END_REF][START_REF] Bayne | Aspects of reproduction in bivalve molluscs[END_REF][START_REF] Pieters | Tissue composition and reproduction of Mytilus edulis in relation to food availability[END_REF][START_REF] Sunila | Reproduction of Mytilus edulis L. (Bivalvia) in a Brackish Water Area, the Gulf of Finland[END_REF][START_REF] Kautsky | Quantitative studies on gonad cycle, fecundity, reproductive output and recruitment in a Baltic Mytilus edulis population[END_REF][START_REF] Lowe | Aspects of gametogenesis in the marine mussel Mytilus edulis[END_REF][START_REF] Sprung | Reproduction and fecundity of the mussel Mytilus edulis at Helgoland (North Sea)[END_REF][START_REF] Villalba | Gametogenic cycle of cultured mussel, Mytilus galloprovincialis, in the bays of Galicia (N.W. Spain)[END_REF][START_REF] Pulfrich | Seasonal variation in the occurrence of planktic bivalve larvae in the Schleswig-Holstein Wadden Sea[END_REF][START_REF] Cáceres-Martínez | Distribution and abundance of mussel (Mytilus galloprovincialis Lmk) larvae and post-larvae in the Ria de Vigo (NW Spain)[END_REF][START_REF] Frantzen | Recruitment of blue mussels, Mytilus edulis L., on suspended collectors in Finnmark, North Norway (70-71°N)[END_REF][START_REF] Duinker | Gonad development and spawning in one and two year old mussels (Mytilus edulis) from Western Norway[END_REF][START_REF] Dias | Gametogenic asynchrony of mussels Mytilus in a mixedspecies area: Implications for management[END_REF][START_REF] Doherty | Synchronous reproduction may facilitate introgression in a hybrid mussel (Mytilus) population[END_REF][START_REF] Troost | Effects of an increasing filter feeder stock on larval abundance in the Oosterschelde estuary (SW Netherlands)[END_REF].

Based on our results, larval presences of mussels appear not to be determined by small-scale (local) variation in environmental conditions (null hypothesis), but by environmental conditions that operate at larger (100s to 1000s of km) spatial scales. The best model (Fig. 6) suggested that, for the North Sea and the Iberian Coast, the probability of occurrence of mussel larvae is high at the end of March with a second peak in end September (North Sea) to the beginning of November (Iberian Coastal). The model results for the Celtic-Biscay and Mediterranean areas were different for the beginning and the end of the year (Fig. 6), which suggested that conditions in the study year and the following year were different. To exclude the possibility of an artefact in the models, this type of study ideally needs to be conducted across more than a single year. The second best but comparably good model (Fig. 7) suggested that on seasonality in occurrence of mussel larvae showed a common seasonal pattern with one peak around mid-April for all sites.

Previous observations on the reproductive cycle of M. edulis suggested a shift from one to two periods of reproduction within a year during the past decades, i.e. the observations of a single spawning event in spring were made in the 1950s to the 1980s [START_REF] Chipperfield | Observations on the breeding and settlement of Mytilus edulis (L.) in British waters[END_REF][START_REF] Pieters | Tissue composition and reproduction of Mytilus edulis in relation to food availability[END_REF][START_REF] Sunila | Reproduction of Mytilus edulis L. (Bivalvia) in a Brackish Water Area, the Gulf of Finland[END_REF][START_REF] Kautsky | Quantitative studies on gonad cycle, fecundity, reproductive output and recruitment in a Baltic Mytilus edulis population[END_REF][START_REF] Sprung | Reproduction and fecundity of the mussel Mytilus edulis at Helgoland (North Sea)[END_REF], whilst multiple spawning in both spring and autumn was observed from the 1990s onwards [START_REF] Villalba | Gametogenic cycle of cultured mussel, Mytilus galloprovincialis, in the bays of Galicia (N.W. Spain)[END_REF][START_REF] Pulfrich | Seasonal variation in the occurrence of planktic bivalve larvae in the Schleswig-Holstein Wadden Sea[END_REF]Cáceras-Martinez and Figueras, 1998;[START_REF] Frantzen | Recruitment of blue mussels, Mytilus edulis L., on suspended collectors in Finnmark, North Norway (70-71°N)[END_REF][START_REF] Duinker | Gonad development and spawning in one and two year old mussels (Mytilus edulis) from Western Norway[END_REF][START_REF] Dias | Gametogenic asynchrony of mussels Mytilus in a mixedspecies area: Implications for management[END_REF][START_REF] Troost | Effects of an increasing filter feeder stock on larval abundance in the Oosterschelde estuary (SW Netherlands)[END_REF]. It is not possible, however, to infer from this information whether the reproductive behaviour of mussels has actually changed during the past decades, because most observations in the older publications were restricted to a limited period of the year.

This possibility of long-term changes in phenology of mussels deserves future investigation, e.g. by studying the phenology of spawning of mussels from various origins under different temperature and food conditions, because climate-induced changes in spawning frequencies have been observed for other marine invertebrates in Europe during the past decades. In south-west Britain, for example, Patella depressa (a warm temperate gastropod species) originally showed single spawning events only in the 1940s whilst multiple periods of spawning and gonad redevelopment occurred during the reproductive season throughout the 2000s [START_REF] Moore | Phenological changes in intertidal con-specific gastropods in response to climate warming[END_REF].

Conclusion

Based on the variety in environmental conditions such as temperature and chlorophyll-a concentrations at our sampling stations (Fig. 2) and potential patchiness in larval distributions, we would have expected that the null hypothesis (local variation) would best explain observed occurrence of oyster and mussel larvae throughout Europe. Both bivalve species, however, exhibited synchronised patterns in occurrence of larvae at mesoscale (100s of km) to synoptic scales (1000s of km). Such synchrony in larval occurrence suggests that environmental conditions relevant to bivalve reproduction and larval survival are more or less similar at these spatial scales.

As the result of our comparative approach, we cannot determine the actual underlying mechanisms for these patterns, nor exclude the possibility that different locations harboured different subpopulations adapted to local conditions. Our results might help, however, to design experiments to test competing hypotheses regarding causal mechanisms such as food availability and temperature at the appropriate ranges and scales [START_REF] Underwood | Structure of intertidal assemblages in different locations: how can local processes be compared[END_REF][START_REF] Jenkins | Spatial and temporal variation in settlement and recruitment of the intertidal barnacle Semibalanus balanoides (L.) (Crustacea: Cirripedia) over a European scale[END_REF]. The unravelling of the nature and strength of structuring factors of phenology of bivalve reproduction is of particular interest in the light of changing environmental conditions as the result of global climate change and the possible consequences for marine food webs and ecosystem services (e.g., [START_REF] Philippart | Climate-related changes in recruitment of the bivalve Macoma balthica[END_REF][START_REF] Kirby | Synergistic effects of climate and fishing in a marine ecosystem[END_REF].

Table 1. Name and positioning of sampling sites, ranked from high to low latitudes, and the name of the Large Marine Ecosystem 1 (www.lme.noaa.gov) and Marine Ecoregion [START_REF] Spalding | Marine ecoregions of the world: A bioregionalization of coastal and shelf areas[END_REF] in which the sites are situated. Furthermore, the potential 2 presence of Mytilus groups is indicated (Mt = M. trossulus, Me = M. edulis, MgA = Atlantic M. galloprovincialis, MgM = 3 Mediterranean M. galloprovincialis) based on literature [START_REF] Quesada | A multilocus allozyme discontinuity in the mussel Mytilus galloprovincialis: the interaction of ecological and life-history factors[END_REF][START_REF] Śmietanka | Mitochondrial DNA lineages in the European populations of mussels (Mytilus spp.)[END_REF], 2009)

is s i is p f a      

H2

There is one underlying latitude effect and one underlying seasonal trend for all the time series logit( )

(DayNr ) (Latitude ) is is is i is p f f a        H3 The latitude effect changes over time logit( ) (DayNr ,Latitude ) is is is i is p f a      

H4

All the time series have four underlying seasonal trends which are related to LMEs

logit( ) (DayNr ) is j is i is p f a       j = 1, ..
( ) (DayNr ) is k is i is p f a       k = 1, ..

5, where k refers to MEs

H6

There is one seasonal trend and the latitude effect differs

per LME logit( ) (DayNr ) (Latitude ) is is j is i is p f f a        j = 1, ..

4, where j refers to LMEs

H7

There is one seasonal trend and the latitude effect differs

per ME logit( ) (DayNr ) (Latitude ) is is k is i is p f f a        k = 1, ..

5, where k refers to MEs

H8

There is one seasonal trend and there is a subspecies effect logit( ) (DayNr )

is is i is p f Subspecies a       

H9

There is one seasonal trend, the latitude effect is the same for all the time series and there is a subspecies effect logit( ) with a 1-dimensional smoother (see Model H1 in Table 4). The vertical axis

(DayNr ) (Latitude ) is is is i is p f f Subspecies a        
shows the fitted values (logistic scale) and the horizontal axis the day number (with values between 1 and 368).

  Crassostrea gigas and Mytilus edulis/galloprovincialis were developed to amplify a part of the cytochrome oxidase subunit I (COXI) mitochondrial gene (Table 2). Primers for C. gigas had at least one mismatch at the 3' terminal base with C. angulata, C. sikamea, C. nippona, C. ariakensis. Primers for M. edulis/galloprovincialis matched exactly with most sequences of female individuals of these species, except for some sequences from Nedlands (Australia) and Dichato (Chile).

  4, where j refers to LMEs H5 All the time series have five underlying seasonal trends which are related to MEs logit

FIGURE CAPTIONS FIGURESFigure 1 .

 CAPTIONS1 FIGURE CAPTIONS

Figure 2 .

 2 Figure 2. Occurrence of larvae of Pacific oysters (Crassostrea gigas),

Figure 3 .Figure 4 .

 34 Figure 3. Estimated smoother obtained by the model in H1 (see Table4) for

Figure 6 .Figure 7 .

 67 Figure 6. Estimated smoothers obtained by the model in H4 for the

  We tested if the observed seasonal patterns of presence of bivalve larvae at 12 European stations could be best explained:

		by local variation dominating the geographical variation in larval
		presence (null hypothesis), or
		by one common seasonality pattern for all stations within Europe
		(Hypothesis 1), or
		by groups of seasonal patterns which are related to specific spatial
		scales such as latitude, large marine ecosystems, and marine
		ecoregions (Hypotheses 2 to 5), or
		by taxonomic groups (Hypothesis 8), or
		by several combinations of these factors (Hypotheses 6, 7 and 9)
		(Table 3).

European marine waters can be classified into several Large Marine Ecosystems, LMEs (Table

1

; Fig.

1

). These LMEs are relatively large regions on the order of 200,000 km 2 or greater, characterized by distinct bathymetry, hydrography, productivity, and trophically dependent populations (see http://www.lme.noaa.gov). Based on a more detailed distribution overview of marine species than the one that was originally used for the definition of the LMEs,

[START_REF] Spalding | Marine ecoregions of the world: A bioregionalization of coastal and shelf areas[END_REF]

subdivided the coastal and shelf areas of the world into 232 marine ecoregions (MEs). Within European marine waters, however, most MEs more or less overlap with LMEs with exception of the Iceland Shelf and the Norwegian Sea (which are each split into 2 ecoregions)

Table 2 .

 2 Primer sequences and additional information 8

	.

Table 3 .

 3 Hypotheses and statistical models for the presence/absence data of larvae of Crassotrea gigas and Mytilus 11 edulis/galloprovincialis in coastal European waters at a synoptic scale. (LME = Large Marine Ecosystem, ME = Marine Ecoregion).

	12
	13
	#	Hypothesis

Table 4 .

 4 Values and differences of the Akaike Information Criterium (AIC), and the Akaike weight (w i ) for models 0 to 9 in Table4. Model. Models for which the level of empirical support for model i is substantial (i.e. the value of Δ i is between 0 and 2) are underlined, the best fit (for which AIC has the lowest value) is printed in bold.

		Crassostrea gigas		Mytilus
					edulis/galloprovincialis
		AIC	Δ i	w i	AIC	Δ i	w i
	H0	119.9	16.0	0.000	281.0	6.6	0.019
	H1	103.9	0.0	0.576	275.9	1.5	0.246
	H2	105.6	1.7	0.246	279.2.	4.8	0.047
	H3	106.3	2.4	0.173	281.1	6.7	0.018
	H4	116.6	12.7	0.001	274.4	0.0	0.519
	H5	120.6	16.7	0.000	277.7	3.3	0.100
	H6	114.6	10.7	0.003	288.0	13.6	0.001
	H7	118.5	14.6	0.000	291.5	17.1	0.000
	H8	not relevant for this species	281.8	7.4	0.013
	H9	not relevant for this species	279.7	5.3	0.037

Acknowledgements

The authors thank the coordinators and funders of the EU-Network of Excellence "Marine Biodiversity and Ecosystem Functioning" (MarBEF) for financially supporting this research. We acknowledge the SOMLIT French network for the provision of environmental data. We are grateful to Jan Beukema, Mardik Leopold and four anonymous referees for their valuable comments on earlier drafts of this manuscript. We are particularly grateful to John S. Gray, who initiated the sampling efforts near Oslo in 2006 but passed away in October 2007.

Table 5. Crassostrea gigas. The estimates, estimated degrees of freedom (obtained by cross-validation) and approximate significance for intercepts and smoothers of best models for oyster larvae (see Table 4). Day number is referring to the day number of the experiment, ranging from 1 (December 15, 2006) to 368 (December 18, 2007). To fit models 2 and 3, both covariates were scaled to have the same variation. 4). The vertical axis shows the fitted values (logistic scale) and the horizontal two axes the day number of the experiment (with values between 1 and 368) and the latitude (H2-LAT; from 45⁰N to 80⁰N). To fit the two-dimensional smoother, both covariates were scaled to have the same variation, but this did not affect the AIC.

Model