# Modelling the spatial heterogeneity of ecological processes in an intertidal estuarine bay: dynamic interactions between bivalves and phytoplankton

K. Grangeré, Sébastien Lefebvre, C. Bacher, P. Cugier, A. Ménesguen

► To cite this version:

K. Grangeré, Sébastien Lefebvre, C. Bacher, P. Cugier, A. Ménesguen. Modelling the spatial heterogeneity of ecological processes in an intertidal estuarine bay: dynamic interactions between bivalves and phytoplankton. Marine Ecology Progress Series, 2010, 415, pp.141-158. 10.3354/meps08659. hal-00765511

HAL Id: hal-00765511

https://hal.science/hal-00765511

Submitted on 19 Dec 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling the spatial heterogeneity of ecological processes in an intertidal estuarine bay: dynamic interactions between bivalves and phytoplankton

Karine Grangeré<fnref n="1" />,<fnref n="2" />*, Sébastien Lefebvre<fnref n="1" />,<fnref n="3" />, Cédric Bacher<fnref n="2" />, Philippe Cugier<fnref n="2" />, Alain Ménesguen<fnref n="2" />

**ABSTRACT:** Spatial patterns in ecological communities result from a combination of physical and biological factors. In an estuarine intertidal bay, spatial differences have been found in the structure of phytoplanktonic communities and in the biological performance of cultivated oysters. It has been hypothesised that trophic heterogeneity exists, although the mechanisms controlling it remain undefined. Spatial and temporal interactions in the structure of phytoplanktonic biomass and in the biological performance of cultivated oysters were highlighted in this estuarine intertidal bay using a 2-dimensional hydrodynamic model coupled to a nutrient-phytoplankton-zooplankton bivalve food web model. The coupled models allowed a reproduction in space and time of variations in the main variables (i.e. nutrients, chlorophyll a (chl a) and bivalve growth and reproduction). Spatial patterns of chl a at the bay scale showed a dichotomy between the eastern and western parts of the bay, with a sharp drop in concentrations above the oyster area. At the smallest scale, significant spatial heterogeneity was obtained in terms of oyster dry weight (DW), with a difference of around 3.0 g between the lowest and the highest oyster DW. Influences of physical and biological factors were discriminated for spatial patterns of phytoplankton at a large scale and for spatial patterns of bivalves at a small scale. Bivalve density, immersion time (i.e. feeding time) and current velocity were identified as the main factors controlling the spatial patterns of phytoplankton and bivalve growth. The results of the model indicate that the effects of spatial scales are much larger than those of temporal scales; this conclusion differs from that expected through observations only. Top-down effects of oysters on phytoplankton biomass at local scales were revealed, whereas bottom-up effects drove primary productivity at the whole bay scale. In general, we conclude that spatial modelling is particularly appropriate to reveal spatial properties which would be difficult to observe directly. Knowledge of ecosystem functioning would be enhanced accordingly.

**KEY WORDS:** Spatial heterogeneity · Hydrodynamics · Ecosystem model · Food supply · Bivalves · Dynamic energy budget model · Physiological status · Baie des Veys

Resale or republication not permitted without written consent of the publisher

## INTRODUCTION

Processes of spatial heterogeneity have been shown to exert a crucial influence on ecosystem function (e.g. Dutilleul 1993, Legendre 1993, Borcard et al. 2004).

Both abiotic and biotic variables that govern the structuring of ecosystems display spatial patterns (Borcard et al. 2004), so that living beings in nature are neither distributed uniformly nor completely at random (Legendre 1993). Indeed, the spatial heterogeneity of a

<footnote n="1">¹Université de Caen Basse-Normandie, UMR 100 IFREMER-UCBN PE2M Physiologie et Ecophysiologie des Mollusques Marins, Esplanade de la Paix, 14032 Caen cedex, France</footnote>

<footnote n="2">²IFREMER, Département Dynamiques de l'Environnement côtier, Technopôle Brest-Iroise, z.i. Pointe du diable, B.P. 70, 29280 Plouzané, France</footnote>

<footnote n="3">³Université de Lille1 sciences et technologies-CNRS UMR 8187 LOG «Laboratoire d'Océanographie et de Géoscience», Station Marine de Wimereux, 28 avenue Foch, 62930 Wimereux, France</footnote>
biological variable may be due not only to underlying environmental constraints (physical or biological), themselves spatially structured by their own generating processes, but also to dynamic processes (e.g. growth, mortality, competition) within these biological variables (Legendre et al. 2002). These spatial structures are usually divided into different classes, such as patches and gradients (Legendre 1993).

Interactions between macrobenthos and their physical and biogeochemical environment have been extensively studied in order to understand the spatial distribution of natural macrobenthic assemblages (e.g. Millet & Guelorget 1994, Legendre et al. 1997, Sousa et al. 2006). The main physical factors influencing the distribution of macrobenthos include salinity, tidal elevation (Fujii 2007), wind-wave disturbance, immersion time (Legendre et al. 1997, Ysebaert et al. 2003), sediment features for deposit feeders, or tidal current velocity for suspension feeders (Herman et al. 1999). The main biotic factor influencing the spatial pattern of benthic biomass is food availability (Fujii 2007, van der Wal et al. 2008), and it is well known that abiotic processes (i.e. hydrodynamics) play a key role in controlling the food supply to benthic filter feeders (Legendre et al. 1997, Ysebaert et al. 2003, Simpson et al. 2007). Conversely, the presence of suspension feeders in an ecosystem often leads to depletion of phytoplankton biomass in the water column (Cloern 1982, Officer et al. 1982, Alpine & Cloern 1992, Prins et al. 1997), which can create significant horizontal and vertical gradients (Simpson et al. 2007). In areas with high benthic biomass, phytoplankton dynamics are strongly linked to benthic processes (Officer et al. 1982, Herman et al. 1999). In addition to these biotic processes, plankton are also structured by abiotic processes such as currents, winds, or vertical stratification/mixing, which drive biological response mechanisms (Hayward & McGowan 1985, Mackas et al. 1985, Ludovisi et al. 2005) and cause the appearance of gradients or patchy structures (Legendre 1993). This is especially the case for continental shelf systems, where the spatial patterns of plankton are usually larger than in the open ocean (Mackas et al. 1985). The study of spatial patterns is thus a crucial step to understanding ecosystem functioning (Borcard et al. 2004).

In the literature, numerous statistical models are commonly used to explore the link between the spatial structuring of biological variables and the structure of the biotic and abiotic environment (e.g. Legendre 1993, Belgrano et al. 1995, Keitt et al. 2002, Legendre et al. 2002, Borcard et al. 2004). These models are based on field surveys and aim to relate biological response variables (e.g. the growth of individuals, the abundance of a species) to explanatory environmental variables (e.g. bathymetry, wind speed, current velocity, predator abundance) (Legendre et al. 2002). However, studies

based on field sampling are subject to some limitations. In order to assess the whole spatial variability it is essential to define an adequate sampling design (Legendre et al. 1997) with sufficient samples, which can be expensive and difficult to implement. Moreover, one of the main limiting factors is the large number of environmental variables that can influence the response mechanisms of the studied biological variable (Legendre et al. 1997). In this context, the spatial patterns studied can be controlled by unmonitored environmental variables (Belgrano et al. 1995), especially in pelagic environments, where spatial structures are weak and transitory (Ludovisi et al. 2005). Approaches based only on statistical models do not allow the dynamics of the spatial pattern of biological variables (i.e. plankton and benthic variables) to be explored, but can provide initial information for the development of numerical models, which are ideally suited to explore the spatio-temporal dynamics of key variables.

Most numerical models in the literature were developed to study spatio-temporal patterns at a large scale (e.g. Moll & Radach 2003 and references therein, Ménesguen et al. 2007), and mainly focus on pelagic spatial patterns (i.e. phyto- and zooplankton, fish). Conversely, only a few numerical modelling studies explicitly take into account the dynamics of both planktonic and benthic species in coastal or estuarine ecosystems (e.g. Bacher et al. 2003, Grant et al. 2008, Duarte et al. 2008), and these models are mainly seston depletion models (Duarte et al. 2008). Bacher et al. (2003) used the coupling between a depletion model, a scope for growth model for *Chlamys farreri* and a 1-dimensional horizontal transport equation. This model could be used to assess the influence of the spatial variation of certain environmental variables (e.g. food availability, temperature, hydrodynamics) on scallop growth. More complex coupled physical-biogeochemical-growth models were mainly developed in order to estimate carrying capacity (Duarte et al. 2003) or to improve the management of bivalve farming (Pastres et al. 2001, Marinov et al. 2007, Spillman et al. 2008, Maar et al. 2009). These models aimed to simulate bivalve growth and did not focus on energy allocation. To the best of our knowledge no previous modelling studies have, therefore, yet been conducted on the physiological response of benthic species facing spatio-temporal variations in biotic and abiotic conditions.

The aim of the present work was to emphasise the eco-physiological response of suspension-feeding bivalves in relation to the spatial heterogeneity of environmental factors and to identify the main physical (e.g. immersion time, hydrodynamics) and biological processes (i.e. food availability) influencing growth patterns and the physiological status of these organisms. Cultivated oysters were chosen as a biological model, because these cultured
populations have the advantage of being regularly mon-
itored, providing numerous spatio-temporal data sets to
calibrate and validate a numerical model. For the same
reason, the Baie des Veys ecosystem (eastern English
Channel, France) was chosen as a study site, because it
is an estuarine bay highly influenced by tides where hy-
drodynamics must have an influence on the spatial dis-
tribution of pelagic variables and on the food availability
for benthic organisms. Furthermore, a previous study
showed that oyster growth in this ecosystem is not homo-
geneous, but exhibits a spatially structured pattern (Cos-
til et al. 2005). However, although this previous study
based on field sampling revealed some spatio-temporal
patterns, it could not explain the entirety of spatial pat-
terns in the bay. In the same way, some recurrent abnor-
mal events of summer oyster mortalities occurred in this
ecosystem. A strong spatial heterogeneity of oyster mor-
talities was recorded at the farm scale (Samain & Mc-
Combie 2008). Thus, the influence of the spatial variabil-
ity of environmental factors on oyster physiology still
needed to be studied.

In the present study, questions were approached on
2 different spatial scales: on the scale of the whole bay
for the spatial pattern of phytoplankton and on the
scale of the oyster-farming area for the spatial pattern
of oyster growth. To this end, a hydrodynamic sub-
model was coupled with an ecosystem sub-model. The
hydrodynamic sub-model was developed with more
detail inside the bay with regards to the oyster-farming
area in order to fully reproduce the spatial pattern at
this smaller scale. The ecosystem sub-model was itself
a coupling between a biogeochemical sub-model,
which simulated the trophic resources of oysters (i.e.
phytoplankton), and an oyster ecophysiological sub-
model, which simulated energy allocation between

growth and reproduction on a mechanis-
tic basis. Results obtained enabled us to
study the magnitude of spatial variabil-
ity in oyster growth and to identify the
factors underlying it. Furthermore, the
results of spatial variability were com-
pared with those of the inter-annual
variability that were extensively dis-
cussed in Grangeré et al. (2009a).
Finally, the definition of an indicator of
energy needs allowed to us to examine
the spatial patterns in oyster physiologi-
cal status.

MATERIALS AND METHODS

**Study area.** Baie des Veys is located
on the French coast of the English
Channel, in the western part of the

Baie de Seine (Fig. 1). It is an estuarine bay under a
macrotidal regime (the maximum tidal range reaches
8 m) with an intertidal zone of around 37 km²
(Desprez et al. 1986). The mean depth is ~4 m and
increases gradually outside of the bay, reaching 25 m
at the sea boundary. In Baie des Veys, the main culti-
vated species is the Pacific oyster *Crassostrea gigas*.
The intertidal area supports considerable oyster-
farming activity (10 500 t of biomass) in an area of
~160 ha located in the north-eastern part of the bay
(Fig. 1). Around 600 t of mussels *Mytilus edulis* are
cultivated in the southern oyster-farming area. Oys-
ters and mussels are also cultivated in the north-
western part of the bay, with standing stocks of
around 3100 and 1330 t, respectively. Wild cockles
*Cerastoderma edule* are found in the southern and
the western parts of the bay in quantities estimated at
1000 and 1500 t, respectively. Freshwater enters the
southern part of the bay from 4 rivers that drain a
watershed area of ~3465 km² (Lefebvre et al. 2009).
The total mean discharge is ~53 m³ s⁻¹ (ranging from
2.2 m³ s⁻¹ at the lowest water level to 400 m³ s⁻¹ dur-
ing floods). The Vire is the main river flowing into the
Baie des Veys and supplies 40% of total input. The
primary production dynamics of the Baie des Veys
ecosystem have been described by Jouenne et al.
(2007). These results indicate that riverine nutrient
inputs lead to a phytoplanktonic bloom during spring,
with diatoms predominating throughout the year.
This bloom is initially limited by a decrease in silicon
availability in late spring, followed by a decrease in
nitrogen concentration ~2 wk later. Generally, *Phaeo-
cystis* cells were seen to develop between the periods
of silicon and nitrogen limitation in spring (Jouenne
et al. 2007).

<figure><img src="image_2.png" /></figure>

Fig. 1. Baie des Veys estuary (modified after Dubois et al. 2007). The 4 points represent the sampling stations for which biogeochemical measurements were available (chl a, nutrients, salinity and temperature). Stns E & S were located near the oyster-farming area in the eastern part of the bay. A third station was located in the north (N) and another in the west (W)
**Model characteristics and implementation in Baie des Veys.** The SiAM 3D model was developed by Cugier & Le Hir (2002). It is a 3-dimensional hydrodynamic model that solves the Navier-Stokes equations by applying the Boussinesq approximation, hydrostatic equilibrium and incompressibility. It is based on a finite difference 'C' grid as defined by Arakawa & Lamb (1977). Since no thermal stratification has been shown in this ecosystem (Jouenne et al. 2007), the model developed for the Baie des Veys was 2-dimensional (2D) and vertically integrated. Haline stratification occurs on some occasions during winter due to high river discharges. However, it did not seem essential to explicitly simulate winter haline stratification considering the aim of the present work. The bottom topography was based on data acquired by the Service Hydrographique et Océanographique de la Marine (SHOM) for the area located outside of the bay and on data from the 'Parc Régional des Marais du Cotentin et du Bessin' for within the bay. The model covers the area between 49°20'N and 49°32'N and between 0°5'W and 1°16'W. The model grid is an irregular rectangular horizontal grid with 11310 squares (130 columns × 87 rows). The spatial resolution was 100 m for the southernmost squares, allowing finer definition locally in the bay and more precisely in the oyster-farming area (i.e. 207 squares) (Fig. 2), then, the mesh size increased regularly towards the northern sea boundaries.

The evolution of both conservative and non-conservative variables was described using an advective-diffusive equation (Eq. 1), where biogeochemical processes provide the terms 'Sources' and 'Sinks':

$$ \frac{\partial(C)}{\partial t} + \frac{\partial(uC)}{\partial x} + \frac{\partial(vC)}{\partial y} = \frac{\partial}{\partial x} \left( k_x \frac{\partial C}{\partial x} \right) + \frac{\partial}{\partial y} \left( k_y \frac{\partial C}{\partial y} \right) + \text{Sources} - \text{Sinks} \quad (1) $$

where C corresponds to any dissolved or particulate variables, u and v represent current velocities in x and y directions (m s⁻¹) and kx and ky correspond to the horizontal diffusion in x and y directions, respectively (m² s⁻¹).

A single sediment layer was taken into account. Bottom exchanges were introduced using a Partheniades formulation for erosion and a Krone formulation for deposition (Cugier & Le Hir 2000). Consequently, the level of suspended particulate matter (SPM) in the water column was the result of river inputs, as well as erosion-deposition processes (Cugier et al. 2005). An adequate simulation of SPM is essential to accurately estimate the light extinction coefficient in these turbid coastal waters. When a grid square has emerged, the calculation of every state variable is stopped until the next immersion.

The model has a land boundary and an ocean one. The forcing variables used were the tidal harmonic components and concentrations for the simulated state variables at the sea boundary, and the light intensity, air temperature, wind speed and cloud cover at the sea surface. Because of the limited available data close to the marine boundary of the model area, state variable values were taken from outputs of the 3-dimensional model of the Baie de Seine developed by Cugier et al. (2005). Tidal harmonic components were provided by the SHOM. Meteorological forcing was provided by

<figure><img src="image_3.png" /><figcaption>Fig. 2. Grid of the hydrodynamic model. (a) The model has irregular grid square sizes allowing different levels of resolution within the bay, with more precision in the oyster-farming area (i.e. squares = 100 × 100 m, inset). (b) The model grid at the scale of the oyster-farming area. P1 to P5: stations where oyster growth was measured</figcaption></figure>
Météo-France at the Englesqueville-lapercée station (north eastern bay). Daily mean measurements (sunshine duration, air temperature, air moisture, atmospheric pressure and cloud cover) or 3 h data (wind speed and direction) were used for the calculation of sea temperature and surface wind-induced stress. Instantaneous solar irradiance was calculated from astronomic calculations corrected for cloudiness (Ménesguen & Hoch 1997). The horizontal grid allows a realistic representation of rivers, respecting channel depths and cross-sections (Cugier & Le Hir 2002, Cugier et al. 2005). River forcings were derived from measurements carried out in the 4 tributaries, from river discharges measured every day by the Direction Régional de l'Environnement (DIREN). Dissolved nutrients (inorganic nitrogen and silicon), suspended particulate inorganic matter (SPiM), Kjeldahl nitrogen and temperature were taken from monthly measurements made by the 'Agence de l'eau Seine-Normandie'. Detrital organic nitrogen concentrations were estimated as Kjeldahl nitrogen minus dissolved inorganic nitrogen (Ménesguen & Hoch 1997). Initial conditions for pelagic variables were determined using field measurements. For benthic variables, initial conditions for nutrients and detrital matter were fixed equal to zero because of the lack of available data.

**The ecosystem model.** The ecosystem model was built from the coupling between a biogeochemical sub-model that simulates the trophic resources of oysters (i.e. primary production and biomass of phytoplankton) depending on environmental conditions (i.e. river inputs, meteorology) and an ecophysiological sub-model that simulates oyster growth and reproduction as functions of trophic resources. This coupling has been fully described by Grangeré et al. (2009a) and applied in a single box model. The following describes the main features of both sub-models.

The biogeochemical sub-model is based on the NPZD model (nutrients, phytoplankton, zooplankton and detritus) developed by Ménesguen & Hoch (1997). The conceptual structure of the ecosystem model is presented in Fig. 3, and state variables are given in Table 1. All differential equations, processes and parameters of nutrients, phyto-

plankton and zooplankton are the same as those calibrated for the Baie de Seine, as fully detailed by Guilaud et al. (2000) and Cugier et al. (2005), with the exception of the primary production equation. Work focussed on the calibration of photosynthetic processes using field measurements of primary production has been performed in separate study (Grangeré et al. 2009b).

A previous study revealed the lack of influence of phosphorus on primary production in the Baie des Veys (Jouenne et al. 2007). In this ecosystem, phytoplanktonic growth is only limited by nitrogen and silicon availability. Thus, the model simulates only the nitrogen and silicon biogeochemical cycles in the

<figure><img src="image_4.png" /><figcaption>Fig. 3. Conceptual diagram of the coupling of the biogeochemical and ecophysiological sub-models for the Baie des Veys estuary. Parameters defined in Table 1</figcaption></figure>

<table><caption>Table 1. State variables of the biogeochemical and ecophysiological sub-models</caption><thead><tr><th>Symbol</th><th>State variable</th><th>Unit</th></tr></thead><tbody><tr><td colspan="3"><strong>Biogeochemical model</strong></td></tr><tr><td>DIN<sub>w</sub></td><td>Pelagic dissolved inorganic nitrogen</td><td>µmol N l<sup>-1</sup></td></tr><tr><td>detN<sub>w</sub></td><td>Pelagic detrital organic nitrogen</td><td>µmol N l<sup>-1</sup></td></tr><tr><td>DIN<sub>sed</sub></td><td>Benthic dissolved inorganic nitrogen</td><td>µmol N l<sup>-1</sup></td></tr><tr><td>detN<sub>sed</sub></td><td>Benthic detrital organic nitrogen</td><td>µmol N l<sup>-1</sup></td></tr><tr><td>Si<sub>w</sub></td><td>Pelagic dissolved inorganic silicon</td><td>µmol Si l<sup>-1</sup></td></tr><tr><td>detSi<sub>w</sub></td><td>Pelagic biogenic detrital silicon</td><td>µmol Si l<sup>-1</sup></td></tr><tr><td>Si<sub>sed</sub></td><td>Benthic dissolved inorganic silicon</td><td>µmol Si l<sup>-1</sup></td></tr><tr><td>detSi<sub>sed</sub></td><td>Benthic biogenic detrital silicon</td><td>µmol Si l<sup>-1</sup></td></tr><tr><td>Diatoms</td><td>Diatoms</td><td>µmol N l<sup>-1</sup></td></tr><tr><td>Zoo</td><td>Zooplankton</td><td>µmol N l<sup>-1</sup></td></tr><tr><td colspan="3"><strong>Ecophysiological model</strong></td></tr><tr><td>E</td><td>Reserve</td><td>J ind.<sup>-1</sup></td></tr><tr><td>V</td><td>Structural body volume</td><td>J ind.<sup>-1</sup></td></tr><tr><td>E<sub>R</sub></td><td>Reproduction buffer</td><td>J ind.<sup>-1</sup></td></tr></tbody></table>
water column and the sediment (Fig. 3). The nitrogen cycle is simulated using 3 state variables: ammonium, nitrate and detrital organic nitrogen. The silicon cycle requires 2 state variables: dissolved inorganic silicon (silicates) and biogenic silicon. The dominance of diatoms in this ecosystem (Jouenne et al. 2007) led us to select only this phytoplanktonic group. Both grazing of diatoms by zooplankton in the water column and consumption of diatoms by wild and cultivated filter-feeders (i.e. cockles, mussels and oysters) were modelled.

The ecophysiological sub-model used was based on dynamic energy budget (DEB) theory (Kooijman 2000). This theory describes the way in which energy derived from food is used for maintenance, growth, development and reproduction (Kooijman 2000). Three state variables are simulated: structural body volume (V), reserves (E) and reproduction buffer (E_R) (Fig. 3, Table 1). All are expressed in energy units (Joules). From these 3 state variables, an additional variable is calculated: dry flesh body weight (g).

Energy assimilated from food is initially integrated in a reserve pool. This energy is a function of the assimilation surface area of the organism and the external food concentration. A fixed proportion of the reserve energy content is spent on growth and maintenance, considering that maintenance has priority over growth. The remaining energy is spent on development in juveniles and on reproduction in adults. The energy allocated to reproduction is stored in the buffer, converted to eggs at the time of reproduction, and emptied at spawning. Spawning occurs when the gonadosomatic index (i.e. the ratio between energy allocated to the reproduction buffer and total energy in the animal) and water temperature reached threshold values. These values are determined as a function of the reproductive cycle of the species concerned. In the present study the ecophysiological sub-model was only applied to oysters, because no available data exist to calibrate and validate mussel or cockle growth.

Growth and reproduction of the Pacific oyster were modelled using the oyster DEB model developed by Pouvreau et al. (2006). In the present study, phytoplankton is the only trophic resource for oysters. In the Baie des Veys ecosystem several sources of organic matter (i.e. phytoplankton, terrestrial and sediment sources) possibly contribute to the oyster diet (Marin-Leal et al. 2008). However, among these different sources, phytoplankton is the main contributor to oyster diet in the Baie des Veys (~50% on a yearly basis). As shown in Grangeré et al. (2009a), the use of phytoplankton alone is sufficient to reproduce the main pattern of oyster dry weight (DW), and, therefore, other sources were not taken into account in the present study.

The coupling between the sub-models consists first of using the outputs from the biogeochemical model (i.e. chl a concentration) as trophic resources for oysters. When this coupling was previously performed in a single box model (Grangeré et al. 2009a), it showed that oyster DW was better simulated using phytoplankton carbon concentrations instead of chl a concentrations as a quantifier for food. Since no carbon measurements were available, phytoplankton carbon concentrations were estimated from simulated chl a concentrations transformed into carbon concentrations using a variable chl a to carbon ratio (chl a:C). To this end, we used the empirical formulation of the chl a:C ratio given by Cloern et al. (1995). Secondly, some processes of the ecophysiological model (i.e. filtration, excretion and biodeposition) provide feedback on nutrient (i.e. dissolved inorganic nitrogen and detrital organic nitrogen) and chl a dynamics simulated by the biogeochemical model. For all bivalves, the filtration rate (FR, see Eq. 2) is expressed in energy units (J d⁻¹). Because the phytoplankton concentration is expressed in chl a units (µg chl a l⁻¹), it was necessary to use a conversion factor between energy units (J) and chl a units. This conversion factor was obtained using a fixed ratio between energy units and phytoplankton carbon concentration (47.76 J mg C l⁻¹; calculated from Platt & Irwin 1973), and the variable chl a to carbon ratio was estimated using Cloern's formulation as explained above.

The formulation of the influence of oysters on the environment (consumption, biodeposition and excretion) has been fully detailed by Grangeré et al. (2009a). In order to adequately simulate the spatio-temporal dynamics of chl a concentration, it was essential to take into account the consumption of all bivalves (i.e. wild and cultivated) on phytoplankton. For all bivalves, the FR was modelled using the same formulation. FR (J d⁻¹) is proportional to the temperature, food and surface area of the organism (Eq. 2, for more information see Grangeré et al. 2009a):

$$FR = IR_{\max} \frac{X}{X_K} V^{2/3} k(T) \quad (2)$$

where X is the food concentration (i.e. phytoplankton carbon concentration µg C l⁻¹), X_K is the half-saturation coefficient (µg C l⁻¹), IR_max is the surface-area-specific maximum ingestion rate (J cm⁻² d⁻¹) and k(T) is the temperature effect described using the Arrhenius equation (Eq. 3; Kooijman 2000).

$$k(T) = \exp\left(\frac{T_A}{T_{\text{opt}}} - \frac{T_A}{T}\right) \times \left[1 + \exp\left(\frac{T_{AL}}{T} - \frac{T_{AL}}{T_L}\right) + \exp\left(\frac{T_{AH}}{T_H} - \frac{T_{AH}}{T}\right)\right] \quad (3)$$
where *T* is ambient temperature, *T*<sub>A</sub> is the Arrhenius temperature, *T*<sub>opt</sub> is the reference temperature, *T*<sub>L</sub> and *T*<sub>H</sub> are the lower and upper boundaries of the temperature tolerance range, respectively, and *T*<sub>AL</sub> and *T*<sub>AH</sub> are the Arrhenius temperatures for the rate of decrease at each boundary. All temperatures are expressed in degrees Kelvin. Parameter values used for the different bivalve species are indicated in Table 2 and come from literature. Only the half-saturation constant for food uptake was free-calibrated. The same value was used in all the simulations.

**Validation data set.** Between May 2002 and October 2003, chl *a* and physicochemical measurements (i.e. water temperature, dissolved inorganic nitrogen [ammonium and nitrate] and silicates) were taken every 2 wk at 4 stations located in the bay (Fig. 1). Stns E and S were located close to the oyster-farming area, Stn N in the north and Stn W in the western part of the bay. These measurements were collected by Jouenne et al. (2007). Considering spatial variability of biogeochemical variables at the scale of the whole bay, comparisons between simulated and measured variables were done for the S and W stations only. These 2 stations were chosen because, in this ecosystem, the main spatial variations (biomass and primary production) occur between the eastern and western parts of the bay (see 'Results: Spatial variability at the scale of the whole bay'). However, in order to adequately validate the ecosystem model, statistical analyses were performed using all available data at the 4 stations.

The only dataset available to validate filter-feeder
growth and reproduction concerned the oysters
located in the north-eastern part of the bay. In this
area, oyster data were collected from January 2002 to
December 2002 at 5 stations (P1 to P5) located within
the oyster-farming area (Fig. 2). Sampling was carried
out every month, except during the reproductive
period (late summer), during which it was done every
2 wk. On each sampling date, 30 individuals were ran-
domly selected, and the dry flesh weight was then
measured at the laboratory. This is the only available
dataset allowing spatial variations in biogeochemical
variables in the bay and spatial variations in oyster
growth to be examined together. Validation of the bio-

geochemical model was thus performed on years 2002
and 2003, and validation of the ecophysiological model
was performed on 2002.

Regression analyses were performed between simu-
lated and measured chl *a*, nitrate, ammonium and
silicate concentrations, temperature and oyster DW.
For each variable we tested whether the determination
coefficient were significantly different from 0.

**Simulations.** All simulated scenarios of spatial vari-
ability of phytoplankton and oyster growth were per-
formed for the years 2002 and 2003. However, due to
the similarity between years for the spatial patterns in
phytoplankton and oyster growth, only the results of
the year 2002 are presented here.

In order to assess whether oyster density had an
influence on the spatial variability of the main biologi-
cal variables (i.e. phytoplankton and oyster growth),
2 scenarios were used. The first one utilised the fully
coupled model described above and cultivated oyster
stock distributed homogeneously in all grid squares
located in the oyster-farming area. The second one
used the same model, but without taking into account
the influence of oyster stock on the environment (i.e.
phytoplankton consumption, excretion and biodeposi-
tion processes); in this scenario, only 1 oyster was
introduced in all grid squares located in the oyster-
farming area. Throughout the present paper these
2 scenarios are referred as Scenario 1 and Scenario 2,
respectively. The comparison between these scenarios
enabled us to study the influence of oyster stock on the
spatial variability of oyster growth.

The spatial and temporal variability of oyster physio-
logical status was studied using an indicator of energy
needs corresponding to the ratio between energy
demand (i.e. maintenance processes) and energy input
(i.e. assimilation), such that:

$$
E_n = \frac{\text{energy demand}}{\text{energy input}} = \frac{\text{maintenance}}{\text{assimilation}}
$$

We assumed that if energy demand was higher than
energy input, then oysters could be in physiological
distress, because neither the food available nor the
level of reserves would be sufficient to sustain mainte-
nance needs.

<table><caption>Table 2. Parameter values used in the calculation of the filtration rate for the various bivalve species found in the Baie des Veys. Data according to Van der Veer et al. (2006). IR<sub>max</sub>: the surface-area-specific maximum ingestion rate (J cm<sup>-2</sup> d<sup>-1</sup>); T<sub>opt</sub>: the reference temperature; T<sub>A</sub>: the Arrhenius temperature; T<sub>L</sub> and T<sub>H</sub>: the lower and upper boundaries of the temperature tolerance range, respectively; T<sub>AL</sub> and T<sub>AH</sub>: the Arrhenius temperatures for the rate of decrease at each boundary</caption>
  <thead>
    <tr>
      <th>Species</th>
      <th>IR<sub>max</sub></th>
      <th>T<sub>opt</sub></th>
      <th>T<sub>A</sub></th>
      <th>T<sub>L</sub></th>
      <th>T<sub>H</sub></th>
      <th>T<sub>AL</sub></th>
      <th>T<sub>AH</sub></th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td>Crassostrea gigas</td>
      <td>560</td>
      <td>293</td>
      <td>5800</td>
      <td>281</td>
      <td>305</td>
      <td>75000</td>
      <td>30000</td>
    </tr>
    <tr>
      <td>Mytilus edulis</td>
      <td>196.8</td>
      <td>289</td>
      <td>5800</td>
      <td>275</td>
      <td>296</td>
      <td>45430</td>
      <td>31376</td>
    </tr>
    <tr>
      <td>Cerastoderma edule</td>
      <td>91.5</td>
      <td>301</td>
      <td>5800</td>
      <td>278</td>
      <td>306</td>
      <td>51154</td>
      <td>47126</td>
    </tr>
  </tbody>
</table>
# RESULTS

## Validation of the coupled model

The model successfully simulated temperature (Fig. 4a,b, Table 3). Maxima occurred in summer, at ~19°C and 21°C for 2002 and 2003, respectively, and minima in winter at ~5°C. There was overall good agreement between simulated and measured chl a concentrations for both stations and both years considered (Fig. 4c,d, Table 3). The timing of the spring bloom was well reproduced, although in 2002 the model failed to simulate its magnitude. Nitrate and silicate concentrations were also adequately simu-

lated. The timing of the decrease in nitrate and silicate concentrations in spring was particularly well reproduced (Fig. 4e–h, Table 3). However, some discrepancies were found at each station, mainly during winter, where the simulation overestimated nitrate and silicate concentrations. A satisfactory simulation of ammonium concentrations was obtained for 2002. In contrast, the simulation for 2003 underestimated ammonium concentrations for both stations (Fig. 4i,j), mainly in autumn, and this led to a low value for the determination coefficient (Table 3).

The time course of oyster DW was simulated using carbon concentrations as a quantifier for food. The half-saturation coefficient for food uptake was cali-

<figure><img src="image_5.png" /><figcaption>Fig. 4. Validation of the biogeochemical model, (a, b) water temperature, (c, d) chl a, (e, f) nitrate, (g, h) silicate and (i, j) ammonium. Black line: model simulations. o: field measurements. Validation was performed for 2 stations: Stn W located in the western part of the bay (a, c, e, g, i) and Stn S located in the eastern part of the bay (b, d, f, h, j). See Fig. 1 for station locations</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Fig. 4 (continued)</figcaption></figure>

brated to 150 µg C l⁻¹. On the whole, oyster DW was adequately simulated for the 5 stations studied (Fig. 5, Table 3, see Fig. 2 for station locations). However, for each station, simulations showed a slight underestimation of oyster DW before spawning and a slight overestimation after spawning. Indeed, the decrease in DW during autumn and winter was not reproduced by the model. The timing of spawning was well simulated for Stns P2, P3 & P4 (Fig. 5), whereas it was less satisfactory for Stns P1 & P5 (Fig. 5). At P1 & P5 spawning seems to extend over several weeks, which is not reproduced by the ecophysiological model. In contrast, the observed spatial differences in oyster DW were adequately simulated by the model. The highest DW was obtained for Stn P4, and the lowest for Stns P1 & P5.

## Spatial variability at the scale of the whole bay

Spatial patterns of all biogeochemical variables were simulated (see 'Results: Validation of the coupled model') and displayed a west to east gradient across the bay; however, only results on the spatial variability of chl *a* are shown in this section. The spatial variability of phytoplankton at the bay scale was studied by taking into account the influence of oysters on the environment (i.e. Scenario 1).

Temporal variability of chl *a* concentration showed that in 2002 the phytoplanktonic spring bloom began outside of the bay, along the Utah beach coast, around the middle of April (see the supplement at www.int-res.com/journals/suppl/m415p141_supp/). The fate of the bloom showed that phytoplankton appeared in the western part of the bay at the beginning of May. The bloom then extended throughout the whole bay and to the Baie de Seine from the middle of May until the end of June, when it decreased rapidly. In Scenario 1, annual mean spatial distribution of chl *a* for 2002 (Fig. 6a) showed maximum values inside the

<table><caption>Table 3. Regression parameters obtained between measured and simulated state variables of the biogeochemical model (chl *a*, nitrate, ammonium, silicates and temperature) and the oyster dynamic energy budget model (dry weight) using a simple linear regression. Determination coefficients were tested according to the equation $Y = aX + b$</caption><thead><tr><th></th><th>R<sup>2</sup></th><th>p-value</th></tr></thead><tbody><tr><td>Chl <i>a</i></td><td>0.49</td><td>&lt;0.001</td></tr><tr><td>Nitrate</td><td>0.25</td><td>&lt;0.001</td></tr><tr><td>Ammonium</td><td>0.11</td><td>&lt;0.001</td></tr><tr><td>Silicates</td><td>0.67</td><td>&lt;0.001</td></tr><tr><td>Temperature</td><td>0.92</td><td>&lt;0.001</td></tr><tr><td>Dry weight</td><td>0.64</td><td>&lt;0.001</td></tr></tbody></table>
bay, ~4 µg chl a l⁻¹. Then, concentrations decreased regularly outside the bay, down to 1 µg chl a l⁻¹ at the oceanic boundaries. Inside the bay, a dichotomy between eastern and western areas was revealed (Fig. 6a). Maximum concentrations of ~4 µg chl a l⁻¹ were found in the western bay, whereas a sharp drop in chl a concentrations was obtained in the eastern bay to ~1.5 to 2 µg chl a l⁻¹. This pattern was not fully highlighted in the validation section (Fig. 4c,d) because the drop in chl a concentrations occurred mainly just above the oyster-growing area, whereas

the sampling station (Stn S; see Fig. 1) was located just outwards of the oyster farms. Outside the bay, spatial patterns displayed a more homogeneous structure. A map of standard deviations (SD) (Fig. 6b) showed high variability in the western bay (~3 to 5 µg chl a l⁻¹), mainly in the intertidal area, and lower variability in the eastern part (~1 to 2 µg chl a l⁻¹). The lowest SD was obtained in the oyster farm area (~1 µg chl a l⁻¹). In Scenario 2, the feedback effect of oysters on the environment was not taken into account. Results showed a homogeneous spatial pattern of chl a

<figure><img src="image_7.png" /></figure>

Fig. 5. *Crassostrea gigas*. Validation of oyster dry weight simulated with the coupled model at 5 stations located in the oyster-farming area (north-eastern part of the bay). Black line: model simulations. o: field measurements. See Fig. 2 for station locations)
concentration inside the bay, with an average value of ~5 µg chl a l⁻¹ (Fig. 6c), i.e. 1 µg chl a l⁻¹ higher than in the scenario with oysters. The associated map of SD (Fig. 6d) showed high variability throughout the whole bay (~4 to 7 µg chl a l⁻¹) even in the south eastern part where no variability occurred with Scenario 1. In order to highlight when spatial discrepancies occurred (Scenario 1), observations of chl a were averaged spatially. Results were compared with the

simulated chl a averaged for the same locations (Fig. 6e). Time courses of average simulated and measured chl a showed similar trends throughout 2002, the main discrepancies were obtained during spring. The comparison of SD showed that the main spatial variability occurred during spring and autumn. This pattern was not fully simulated by the coupled model; the magnitude of variability was always lower in the simulation.

<figure><img src="image_8.png" /><figcaption>Fig. 6. Spatial and temporal variations in chl a concentrations for the year 2002. (a,c) Annual mean chl a concentrations and (b,d) associated SD for the whole bay in the simulation (a,b) with oysters *Crassostrea gigas* and (c,d) without oysters, respectively. (e) Time course of mean measured chl a and the corresponding mean simulated chl a. Means ± SD correspond to the 4 stations presented in Fig. 1</figcaption></figure>
Spatial variability at the scale of the oyster-growing
area

In order to quantify the importance of spatial differ-
ences in oyster DW over the year, descriptive statistics
were estimated for the whole grid squares where oys-
ters were present (see Fig. 2 for oyster farmed model
grid). To this end, every day we extracted the median
as well as the 5th, the 25th, the 75th and the 95th per-
centile among the 207 squares located on the oyster
farms (Fig. 7). Results showed a marked difference
between the highest and the lowest values. In all
squares, the initial condition was fixed at 0.4 g. The
lowest final DW ~1.2 g, whereas the highest final DW
~4.0 g. However, all values were not distributed regu-
larly between both extremities; 75% of the individual
oyster final DW were between 1.2 and 2.55 g, and only
20% were between 2.55 and 4.0 g.

In Scenario 1 (i.e. feedback of oyster on phytoplank-
ton biomass through filtration and nitrogen excretion;
Fig. 8b), the highest DW were obtained for oysters
located further from the coast, in the northern part of
the oyster-farming area (~4.0 g). Then, the DW regu-
larly decreased from the open water towards the coast,
giving the lowest values along the coast in the south-
ern part of the oyster farming (~1.0 g). In contrast,
Scenario 2 displayed a different spatial structure
(Fig. 8c). The highest DW obtained in the north-west
(~8.0 g), whereas the lowest DW were found in the
south (~4.0 g). In Scenario 2, the spatial pattern of
oyster DW was very close to the bathymetric spatial
pattern (Fig. 8a). Comparison between final DW and
bathymetry showed the existence of a linear relation-

ship between these 2 parameters (Fig. 8e). This was
not the case with Scenario 1, which considered the
feedback effect of oysters on the environment (Fig. 8d).

In Scenario 2, at a given bathymetric level, varia-
bility in oyster DW was obtained. This could reach
around 0.7 g for the highest bathymetric level (Fig. 8e).
In parallel, heterogeneity in mean Eulerian residual
currents was simulated at the scale of the oyster area
(Fig. 9). The highest current velocities (~0.01 m s⁻¹)
were found in the northern and the southern parts of
the oyster-farming area, whereas the lowest current
velocities (~0.002 m s⁻¹) were found along the coast
and in the central part of the oyster-farming area.
In the latter area, the presence of an eddy with low
current velocities was revealed (Fig. 9).

Spatial variability in the indicator of energy needs ($E_n$) was seen on the oyster farms (Fig. 10). In order to compare the modelled energy need and the observed spatial mortalities, spatial variations of $E_n$ were studied for the middle of September period (i.e. the period of highest mortalities in the Baie des Veys; Fig. 10a). Results displayed the highest values of $E_n$ in the south-ern part of the oyster area (~2 to 2.5) and the lowest values in the north (~0.5 to 1). In this ecosystem a map of oyster mortalities was performed at the farm scale (Samain & McCombie 2008). This map of oyster mor-talities showed a patch of higher mortality in the south-ern part of the oyster area and a decreasing gradient of mortality in the north (Fig. 10b). This pattern is re-current in this ecosystem and occurred every year, the only difference between years is the magnitude of mortalities.

<figure><img src="image_9.png" /><figcaption>Fig. 7. *Crassostrea gigas*. Variability of oyster dry weight at the oyster farm scale over the year 2002. All descriptive statistics were estimated for the whole oyster farming area every day. Corresponding statistics are the 5th and the 25th percentile, the median, the 75th and the 95th percentile, respectively</figcaption></figure>

DISCUSSION

Validation of the coupled model

Overall, the coupled model was able to reproduce
the main space and time patterns for the biogeochem-
ical variables and for oyster *Crassostrea gigas* growth.
The main spatial differences in measurements were
observed for chl *a* concentrations. For both years stud-
ied, the chl *a* level was always higher in the western
part of the bay, and simulations adequately repro-
duced this spatial pattern. In 2002, however, the
model under-estimated the magnitude of the spring
bloom for both stations. In this ecosystem, *Phaeocystis*
is commonly observed to develop in late spring,
between the silicon and nitrogen limitations (Jouenne
et al. 2007). The high chl *a* concentrations observed
during spring 2002 can be attributed to the presence
of *Phaeocystis*, which is not simulated by the biogeo-
chemical model. Jouenne et al. (2007) have observed
that the cell density of *Phaeocystis* is 4 times lower in
2003 than in 2002, and our results were in accordance with this observation because the simulation of diatoms for 2003 was sufficient to reproduce the magnitude of the spring bloom. Several studies have reported that the influence of *Phaeocystis* on bivalve physiology depends on the colony concentration (Smaal & Twisk 1997). A decrease of clearance rate is commonly observed with the increase of colony concentrations. Because colony form dominates during a bloom (Schoemann et al. 2005), it does not seem to be essential to explicitly simulate the *Phaeocystis* bloom in order to improve oyster growth, except while tak-

ing into account the negative effect on filtration and ingestion rate. In another study, Grangeré et al. (2009a) applied the same coupled model on a box model of the Baie des Veys. They obtained the same discrepancies between simulated and measured nutrient concentrations. More particularly, the model did not reproduce the increase of ammonium concentrations in late 2003. Thus, the introduction of spatial variability does not improve the simulation of ammonium concentrations. This can be explained by an under-estimation of mineralisation processes in the sediment or in the water column. Marine ecosystem

<figure><img src="image_10.png" /><figcaption>Fig. 8. *Crassostrea gigas*. Spatial pattern of oyster final dry weight (DW) at the oyster-growing area for the year 2002. (a) Bathymetry of the oyster-farming area. Oyster final DW in simulations (b) with and (c) without oyster feedback. Relation between oyster final DW and the bathymetric level for simulations (d) with and (e) without oyster feedback</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Residual Eulerian currents within the oyster area</figcaption></figure>

models, including the microbial loop, are still scarce. In the MIRO ecosystem model, Lancelot et al. (2005) included the microbial loop with different classes of biodegradability; this improves the simulation of organic matter mineralisation in the water column and the sediment. However, the implementation of this compartment strongly increases the complexity of the ecosystem model. Furthermore, validation of the microbial loop requires available *in situ* data, which are not commonly measured.

Simulation of oyster growth was validated for 5 stations located in various places in the cultivation area. The model was capable of reproducing the observed

spatial variations of oyster growth. According to Grangeré et al. (2009a), in this ecosystem the use of phytoplankton as the only trophic resource for oysters is sufficient to reproduce the main patterns of oyster growth. The main discrepancy in growth between these studies concerns spawning. Although, the timing of spawning is well simulated overall, the analysis of oyster growth measurements shows that, for some locations, spawning can extend over several weeks, which was not reproduced by the model. Two processes can explain these results: (1) In this ecosystem, it is well known that partial spawning can occur (Costil et al. 2005, Enríquez-Díaz et al. 2009), which is not included in the ecophysiological model as only total spawning was simulated (Pouvreau et al. 2006). (2) The extension of spawning could be explained

by an inter-individual variability in the oyster reproductive effort. In the present study, we only investigated the response of an average individual to variations in physical and biological processes; thus, we did not take into account inter-individual variability.

## Spatial pattern at the scale of the whole bay

In coastal ecosystems, spatial structure of the phytoplankton biomass is strongly regulated by physical processes (i.e. currents, winds; Mackas et al. 1985) and biological processes (i.e. reproduction and death,

<figure><img src="image_12.png" /><figcaption>Fig. 10. *Crassostrea gigas*. (a) Spatial pattern of the indicator of energy needs ($E_n$) on 15 September 2002. (b) The case of the Baie des Veys: spatial distribution of oyster mortalities (%) (illustrated by the isolines) from 2000 to 2001, data according to Samain & McCombie (2008, p 43); figure printed by permission (Sylvand et al. 2003; their Fig. 1.31B)</figcaption></figure>
predator–prey interactions; Legendre 1993). In the Baie des Veys ecosystem spatial patterns of phytoplankton are the result of a combination of physical and biological processes occurring at 2 different spatial scales. The development and fate of phytoplankton is mainly controlled by physical processes at a large scale (see the supplement at www.int-res.com/journals/suppl/m415p141_supp/) and biological processes (oyster filtration) at a more local scale (i.e. inside the bay). For the latter, the spatial pattern of chl a concentration showed a strong depletion above the cultivation area and all along the east coast of the bay up to the river mouth. Decrease in phytoplankton biomass is commonly attributed to high consumption rates of natural (Cloern 1982, Officer et al. 1982, Hily 1991, Herman et al. 1999) or cultivated (Grant et al. 2008) suspension feeders. In our study, simulations performed without the influence of oysters on the environment (i.e. Scenario 2) showed a homogeneous chl a concentration throughout the whole bay (Fig. 6c) and a level slightly higher than in simulations with oysters. Thus, we concluded that the sharp drop in phytoplankton biomass was a consequence of the presence of cultivated oysters, while the presence of natural benthos in the other part of the bay did not impact the model output as much. Furthermore, the map of the SD in chl a values indicated low variability in the oyster-farming area, whereas high variability was simulated in the western part of the bay. These results emphasise the top-down control of oysters on phytoplankton biomass at a local scale (i.e. above the oyster-growing area). In other parts of the bay, where physical rather than biological (natural benthos) processes dominate, bottom-up control occurs. This latter point reveals some differences in ecosystem functioning between the oyster-farming area and the whole bay. In the absence of oysters the whole bay would be controlled by bottom-up processes (i.e. nutrient concentrations), indicating that this estuarine ecosystem is strongly influenced by freshwater inputs. The presence of a dense population of cultivated bivalves, however, strongly modified the functioning of the ecosystem, locally with some consequences for the whole bay (e.g. a small drop of phytoplankton biomass).

## Spatial pattern at the scale of the oyster-growing area

As for spatial pattern at the bay scale, the spatial pattern within the oyster-farming area also results from a combination of physical and biological processes. In the preceding sub-section, we highlighted that hydrodynamics have an influence on phytoplankton spatial heterogeneity. At the small scale of the culti-

vated area, this means that physical processes support the food availability for oysters (Legendre et al. 1997, Ysebaert et al. 2003). In order to identify the physical factors that influence the spatial pattern of oyster growth, the effect of food depletion was removed (i.e. Scenario 2). In this scenario, spatial patterns of oyster growth were very close to those of the bathymetric level within the oyster-cultivation area. This indicates that immersion time has an influence on the spatial heterogeneity of oyster growth at this local scale (Legendre et al. 1997). Previous experiments showed that long immersion periods, allowing long oyster feeding times, supported higher body mass (Honkoop & Beukema 1997). Moreover, in the Baie des Veys ecosystem, some spatial differences observed in oyster growth were attributed to immersion time (Costil et al. 2005). However, for given immersion levels, significant variability was shown in final oyster DW. Thus, immersion time did not completely explain the heterogeneity of oyster growth entirely. Another physical factor identified was the current velocity (Legendre et al. 1997, Herman et al. 1999). Indeed, residual Eulerian currents did not display homogeneous structure at the scale of the oyster-growing area. For the scallop *Chlamys farreri*, a positive relationship was found between DW and maximum current velocity (Bacher et al. 2003). However, in order to best explore the relation between current velocity and oyster growth it will be necessary to take into account the influence of oyster culture structures (iron racks with their feet buried in the sediment), which could significantly modify the current circulation pattern. In the intertidal area, food availability for bivalves is strongly regulated by hydrodynamic processes and their interaction with bathymetry.

At the scale of the whole bay, a significant depletion of phytoplankton was simulated above the oyster-cultivated area. This indicates that oyster density has an influence on ecosystem functioning. The large difference obtained in the spatial patterns of oyster final DW between Scenarios 1 and 2 led us to conclude that cultivated oyster density itself acts as a negative feedback controlling the spatial heterogeneity of oyster growth. The Baie des Veys ecosystem is dominated by tidal flats in such a way that at low tide, depending on the tidal level, most of the bay may be exposed. Tidal circulation in the oyster area is directed from the north towards the south. This implies that available food in the water column is initially filtered by oysters located in the northern and north-western parts of the oyster-farming area. Consequently, for oysters located in the southern part of the farming area and along the coast, the probability of filtering previously phytoplankton-depleted water is greater if the local biomass is high (Herman et al. 1999). This structure is consistent with the depletion of phytoplankton found above the oyster
area. For natural populations, the depletion of food can result in lower growth rates for the bivalves located in the inner patches of a bed (Prins et al. 1997). In another cultivated species, a decrease in scallop DW was highlighted in relation with the increase of food depletion in Sungo Bay in China (Bacher et al. 2003). In the present study, lower growth rates were not seen for oysters located in the inner part of the oyster-farming area, but were confined to oysters grown along the south coast. For these oysters, lower growth rates are the result of a combination between local tidal circulation and intraspecific competition occurring at high oyster density (Prins et al. 1997, Newell 2004).

Modelling spatial processes is critical

Small-scale spatial heterogeneity must be taken into account in order to fully understand ecosystem functioning, but has rarely been investigated (Dubois et al. 2007). Ecosystem models found in the literature explore the spatial heterogeneity of the environment at a large scale (i.e. Moll & Radach 2003, Ménesguen et al. 2007). In contrast, ecosystem models on spatial variability at a small scale remain scarce. Grangeré et al. (2009a) applied an ecosystem model the same as that in the present study to a single box model of the Baie des Veys in order to study the influence of the inter-annual variability in environmental factors (i.e. river nutrient inputs and meteorology) on the biological performances (i.e. growth and physiological status) of the Pacific oyster. Results showed that the inter-annual variability in food availability and water temperature led to shifts in the growth and reproductive patterns of oysters (i.e. initiation of the spring growth period, timing of spawning). However, inter-annual variability did not exceed 0.4 g in terms of final DW between wet and dry years. In contrast, the magnitude of the spatial variability highlighted in the present study showed a difference of ~3.0 g between the highest and the lowest final oyster DW obtained during the same year. Spatial variability is at least 3 times greater than inter-annual variability. Studies performed without taking into account the spatial heterogeneity of the environment neglect the major part of variability. Also, the inter-annual variability in oyster physiological status was estimated in terms of energy needs in order to emphasise some periods of potential physiological distress. Recurrent events of abnormal oyster summer mortalities have been observed in this ecosystem (Samain & McCombie 2008). Previous observations showed that these events varied on both spatial and temporal scales. The inter-annual variability in mortality events has been extensively discussed by Grangeré et al. (2009a), highlighting a differential response be-

tween wet and dry years. In the present work, the same indicator was applied in order to explore the spatial patterns of oyster physiological status within the oyster area. In September (i.e. the period of highest mortalities in the Baie des Veys), spatial variations of the indicator of energy needs showed that energy demand is higher than energy input for oysters located in the southern part of the cultivation area only. This pattern is consistent with the recorded spatial heterogeneity in oyster mortalities (Samain & McCombie 2008). Although the link between energy needs and mortalities was not direct, this result allows one to assume that physiological distress could be occurring during this period. This stress could make oysters more sensitive to external factors like pollution or pathogens, which could lead to mortality events. To the best of our knowledge, there has been no previous study aimed to describe the spatial heterogeneity of bivalve energy needs. Even if it is a first step, there is a real need to take into account ecophysiological processes in ecosystem models. This indicator could, therefore, be very useful to identify areas of potential physiological distress in relation to spatial variations of environmental factors. As a whole, the comparison between modelling and mortality studies highlighted that in this ecosystem, the spatial variability scale also seems to be broader than the inter-annual one. These results highlight that spatial heterogeneity is an essential component of ecosystems (Legendre 1993). To this end, ecosystem modelling provides a powerful tool to improve the knowledge of ecosystem variability at different spatial and temporal scales.

CONCLUSIONS

The present study was undertaken in order to identify the key processes that best explain the observed spatial heterogeneity in phytoplankton biomass and oyster growth. Hydrodynamics (currents and winds), immersion time and oyster density were identified as the main factors controlling spatial heterogeneity at both large and small scales. Our results showed that, to simultaneously assess inter-annual and spatial variability, the use of a coupled hydrodynamical–biogeochemical–ecophysiological model is necessary. The model developed in this work could also be used in a pilot study to define sampling plans or field experiments. Simulations could be performed in order to determine the environmental factors involved in the spatial pattern of the studied variables or to define the best location for sampling stations in relation to the spatial pattern of structuring variables. We demonstrated that cultivated oysters have a significant influence on ecosystem functioning, which does not seem to
be the case for natural populations. A future study could be made to attempt to identify the ecological consequences that a modification in human activities would have for natural populations.

*Acknowledgements.* This work was supported by IFREMER and by the Conseil Régional de Basse-Normandie (France). The authors thank Dr. Helen McCombie-Boudry for English language assistance and 3 anonymous reviewers for helpful comments on the manuscript.

## LITERATURE CITED

* Alpine AE, Cloern JE (1992) Trophic interactions and direct physical effects control phytoplankton biomass and production in an estuary. Limnol Oceanogr 37:946-955

Arakawa A, Lamb VR (1977) Computational design of the basic dynamical process of the UCLA general circulation model. Methods Comput Phys 17:173-265

* Bacher C, Grant J, Hawkins AJS, Jianguang F, Mingyuan Z, Besnard M (2003) Modelling the effect of food depletion on scallop growth in Sungo Bay (China). Aquat Living Resour 16:10-24

* Belgrano A, Legendre P, Dewarumez JM, Frontier S (1995) Spatial structure and ecological variation of meroplankton on the French-Belgian coast of the North Sea. Mar Ecol Prog Ser 128:43-50

* Borcard D, Legendre P, Avois-Jacquet C, Tuomisto H (2004) Dissecting the spatial structure of ecological data at multiple scales. Ecology 85:1826-1832

* Cloern JE (1982) Does the benthos control phytoplankton biomass in South San Francisco Bay? Mar Ecol Prog Ser 9: 191-202

* Cloern JE, Grenz C, Vidergar-Lucas L (1995) An empirical model of the phytoplankton chlorophyll: carbon ratio-the conversion factor between productivity and growth rate. Limnol Oceanogr 40:1313-1321

* Costil K, Royer J, Ropert M, Soletchnik P, Mathieu M (2005) Spatio-temporal variations in biological performances and summer mortality of the Pacific oyster *Crassostrea gigas* in Normandy (France). Helgol Mar Res 59:286-300

* Cugier P, Le Hir P (2002) Development of a 3D hydrodynamic model for coastal ecosystem modelling. Application to the plume of the Seine River (France). Estuar Coast Shelf Sci 55:673-695

Cugier P, Le Hir P (2000) Modélisation 3D des matières en suspension en Baie de Seine Orientale (Manche, France). CR Acad Sci Paris 331:287-294

* Cugier P, Ménesguen A, Guillaud JF (2005) Three-dimensional (3D) ecological modelling of the Bay of Seine (English Channel, France). J Sea Res 54:104-124

* Desprez M, Ducrotoy JP, Sylvand B (1986) Fluctuations naturelles et évolution artificielle des biocénoses macrozoobenthiques intertidales de trois estuaires des côtes françaises de la Manche. Hydrobiologia 142:249-270

* Duarte P, Meneses R, Hawkins AJS, Zhu M, Fang J, Grant J (2003) Mathematical modelling to assess the carrying capacity for multi-species culture within coastal waters. Ecol Model 168:109-143

* Duarte P, Labarta U, Fernandez-Reiriz MJ (2008) Modelling local food depletion effects in mussel rafts of Galician Rias. Aquaculture 274:300-312

* Dubois S, Orvain F, Marin-Leal JC, Ropert M, Lefebvre S (2007) Small-scale spatial variability of food partitioning between cultivated oysters and associated suspension-

feeding species, as revealed by stable isotopes. Mar Ecol Prog Ser 336:151-160

* Doutilleul P (1993) Spatial heterogeneity and the design of ecological field experiments. Ecology 74:1646-1658

* Enríquez-Díaz M, Pouvreau S, Chavez-Villalba J, Le Pennec M (2009) Gametogenesis, reproductive investment, and spawning behavior of the Pacific giant oyster *Crassostrea gigas*: evidence of an environment-dependent strategy. Aquacult Int 17:491-506

* Fujii T (2007) Spatial patterns of benthic macrofauna in relation to environmental variables in an intertidal habitat in the Humber estuary, UK: developing a tool for estuarine shoreline management. Estuar Coast Shelf Sci 75:101-119

* Grangeré K, Ménesguen A, Lefebvre S, Bacher C, Pouvreau S (2009a) Modelling the influence of environmental factors on the physiological status of the Pacific oyster *Crassostrea gigas* in an estuarine embayment; the Baie des Veys (France). J Sea Res 62:147-158

* Grangeré K, Lefebvre S, Ménesguen A, Jouenne F (2009b) On the interest of using field primary production data to calibrate phytoplankton rate processes in ecosystem models. Estuar Coast Shelf Sci 81:169-178

* Grant J, Bacher C, Cranford PJ, Guyondet T, Carreau M (2008) A spatially explicit ecosystem model of seston depletion in dense mussel culture. J Mar Syst 73:155-168

* Guillaud JF, Andrieux F, Ménesguen A (2000) Biogeochemical modelling in the Bay of Seine (France): an improvement by introducing phosphorus in nutrient cycles. J Mar Syst 25:369-386

* Hayward TL, McGowan JA (1985) Spatial patterns of chlorophyll, primary production, macrozooplankton biomass, and physical structure in the central north Pacific Ocean. J Plankton Res 7:147-167

* Herman PMJ, Middelburg JJ, Van De Koppel J, Heip CHR (1999) Ecology of estuarine macrobenthos. Adv Ecol Res 29:195-240

* Hily C (1991) Is the activity of benthic suspension feeders a factor controlling water quality in the Bay of Brest? Mar Ecol Prog Ser 69:179-188

* Honkoop PJC, Beukema JJ (1997) Loss of body mass in winter in three intertidal bivalve species: an experimental and observational study of the interacting effects between water temperature, feeding time and feeding behaviour. J Exp Mar Biol Ecol 212:277-297

* Jouenne F, Lefebvre S, Veron B, Lagadeuc Y (2007) Phytoplankton community structure and primary production in small intertidal estuarine-bay ecosystem (eastern English Channel, France). Mar Biol 151:805-825

* Keitt TH, Bjornstad ON, Dixon PM, Citron-Pousty S (2002) Accounting for spatial pattern when modeling organism-environment interactions. Ecography 25:616-625

Kooijman SALM (2000) Dynamic energy and mass budgets in biological systems. Cambridge University Press, Cambridge

* Lancelot C, Spitz Y, Gypens N, Ruddick K and others (2005) Modelling diatom and *Phaeocystis* blooms and nutrient cycles in the Southern Bight of the North Sea: the MIRO model. Mar Ecol Prog Ser 289:63-78

* Lefebvre S, Harma C, Blin JL (2009) Trophic typology of coastal ecosystems based on δ¹³C and δ¹⁵N ratios in an opportunistic suspension feeder. Mar Ecol Prog Ser 390:27-37

* Legendre P (1993) Spatial autocorrelation: Trouble or new paradigm? Ecology 74:1659-1673

* Legendre P, Thrush SF, Cummings VJ, Dayton PK and others (1997) Spatial structure of bivalves in a sandflat: scale and generating processes. J Exp Mar Biol Ecol 216:99-128

* Legendre P, Dale MRT, Fortin MJ, Gurevitch J, Hohn M, Myers D (2002) The consequences of spatial structure for
the design and analysis of ecological field surveys. Ecol-
graphy 25:601–615

► Ludovisi A, Minozzo M, Pandolfi P, Taticchi MI (2005) Modeling the horizontal spatial structure of planktonic community in Lake Trasimeno (Umbria, Italy) using multivariate geostatistical methods. *Ecol Model* 181:247–262

► Maar M, Bolding K, Petersen JK, Hansen JLS, Timmermann K (2009) Local effects of blue mussels around turbine foundations in an ecosystem model of Nysted off-shore wind farm, Denmark. *J Sea Res* 62:159–174

Mackas DL, Denman KL, Abbott MR (1985) Plankton patchiness—biology in the physical vernacular. *Bull Mar Sci* 37: 652–674

► Marín-Leal JC, Dubois S, Orvain F, Galois R and others (2008) Stable isotopes (δ¹³C, δ¹⁵N) and modelling as tools to estimate the trophic ecology of cultivated oysters in two contrasting environments. *Mar Biol* 153:673–688

► Marinov D, Galbiati L, Giordani G, Viaroli P, Norro A, Ben-
civelli S, Zaldivar JM (2007) An integrated modelling
approach for the management of clam farming in coastal
lagoons. Aquaculture 269:306–320

► Ménesguen A, Hoch T (1997) Modelling the biogeochemical cycles of elements limiting primary production in the English Channel. 1. Role of thermohaline stratification. Mar Ecol Prog Ser 146:173–188

► Ménesguen A, Cugier P, Loyer S, Vanhoutte-Brunier A, Hoch T, Guillaud JF, Gohin F (2007) Two- or three-layered box-models versus fine 3D models for coastal ecological modelling? A comparative study in the English Channel (western Europe). *J Mar Syst* 64:47–65

► Millet B, Guelorget O (1994) Spatial and seasonal variability in the relationships between benthic communities and physical environment in a lagoon ecosystem. Mar Ecol Prog Ser 108:161–174

► Moll A, Radach G (2003) Review of three-dimensional eco-
logical modelling related to the North Sea shelf system—
Part 1: models and their results. Prog Oceanogr 57:
175–217

Newell RIE (2004) Ecosystem influences of natural and culti-
vated populations of suspension-feeding bivalve molluscs:
a review. *J Shellfish* Res 23:51–61

► Officer CB, Smayda TJ, Mann R (1982) Benthic filter feeding: a natural eutrophication control. Mar Ecol Prog Ser 9: 203–210

► Pastres R, Solidoro C, Cossarini G, Canu DM, Dejak C (2001) Managing the rearing of *Tapes philippinarum* in the lagoon of Venice: a decision support system. Ecol Model 138:231–245

► Platt T, Irwin B (1973) Caloric content of phytoplankton. Limnol Oceanogr 18:306–310

► Pouvreau S, Bourles Y, Lefebvre S, Gangnery A, Alunno-Bruscia M (2006) Application of a dynamic energy budget model to the Pacific oyster, *Crassostrea gigas*, reared under various environmental conditions. *J Sea Res* 56: 156–167

► Prins TC, Smaal AC, Dame RF (1997) A review of the feed-backs between bivalve grazing and ecosystem processes. *Aquat Ecol* 31:349–359

Samain JF, McCombie H (eds) (2008) Summer mortality of Pacific oyster *Crassostrea gigas*. The Morest Project. Collection Synthèses, Éditions Quae/Ifremer, Plouzané Cedex

► Schoemann V, Becquevort S, Stefels J, Rousseau V, Lancelot C (2005) *Phaeocystis* blooms in the global ocean and their controlling mechanisms: a review. *J Sea Res* 53:43–66

► Simpson JH, Berx B, Gascoigne J, Saurel C (2007) The inter-action of tidal advection, diffusion and mussel filtration in a tidal channel. *J Mar Syst* 68:556–568

► Smaal AC, Twisk F (1997) Filtration and absorption of *Phaeo-cystis globosa* by the mussel *Mytilus edulis* L. *J Exp Mar Biol Ecol* 209:33–46

► Sousa R, Dias S, Antunes JC (2006) Spatial subtidal macrobenthic distribution in relation to abiotic conditions in the Lima estuary, NW of Portugal. *Hydrobiologia* 559: 135-148

► Spillman CM, Hamilton DP, Hipsey MR, Imberger J (2008) A spatially resolved model of seasonal variations in phytoplankton and clam (*Tapes philippinarum*) biomass in Barbamarco Lagoon, Italy. *Estuar Coast Shelf Sci* 79: 187-203

Sylvand B, Marion C, Lecouturier A, Ropert M (2003) Nou-
velle cartographie sédimentaire de la zone conchylicole
de Grandcamp-Maisy en Baie des Veys (Baie de Seine
occidentale, Manche orientale). GEMEL. Programme
Morest WP6 – contrat 02 6 522 053

van der Veer HW, Cardoso J, van der Meer J (2006) The esti-
mation of DEB parameters for various Northeast Atlantic
bivalve species. *J Sea Res* 56:107–124

► van der Wal D, Herman PMJ, Forster RM, Ysebaert T and others (2008) Distribution and dynamics of intertidal macrobenthos predicted from remote sensing: response to microphytobenthos and environment. Mar Ecol Prog Ser 367:57–72

► Ysebaert T, Herman PMJ, Meire P, Craeymeersch J, Verbeek H, Heip CHR (2003) Large-scale spatial patterns in estuaries: estuarine macrobenthic communities in the Schelde estuary, NW Europe. Estuar Coast Shelf Sci 57:335–355