# Ecosystem modelling to assess the impact of rearing density, environment variability and mortality on oyster production

P Cugier, Y Thomas, C Bacher

► To cite this version:

P Cugier, Y Thomas, C Bacher. Ecosystem modelling to assess the impact of rearing density, environment variability and mortality on oyster production. Aquaculture Environment Interactions, 2022, 14, pp.53-70. 10.3354/aei00428. hal-03836435

HAL Id: hal-03836435

https://hal.science/hal-03836435

Submitted on 2 Nov 2022

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
<figure><img src="image_2.png" /></figure>

FEATURE ARTICLE

OPEN
ACCESS

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

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

# Ecosystem modelling to assess the impact of rearing density, environment variability and mortality on oyster production

Philippe Cugier<fnref n="1" />*, Yoann Thomas<fnref n="2" />, Cédric Bacher<fnref n="1" />

ABSTRACT: The natural productivity of ecosystems, farming practices and mortality events drive the rearing density and growth of oysters in shellfish farming areas. The variability of these drivers, which can be of natural or anthropic origin, is therefore an important source of variation in the growth performance and production of shellfish. Knowledge of these variables and their relative importance help producers to anticipate their effects and adapt aquaculture practices in order to limit negative impacts and guarantee a constant, or at least acceptable, level of production. In this paper, we implement a 3D model coupling hydrodynamics, primary production and individual growth to predict oyster growth and production in Bourgneuf Bay (French Atlantic coast). We set up numerous scenarios to compare and hierarchize the impacts of aquaculture practices, environmental variability and mortality events on shellfish production. Our results allowed us to propose a simplified management tool, in the form of response functions, to optimize shellfish farming practices interannually. This tool will help shellfish farmers quickly recover production levels in response to variations in mortality and/or environmental conditions.

KEY WORDS: 3D ecosystem modelling · Individual growth model · Oyster · Environment variability · Mortality · Scenarios · Management tool

## 1. INTRODUCTION

With an annual production of around 80 million t, among which 17 million t are represented by molluscs (FAO 2018), aquaculture contributes significantly

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

to the growing food needs of human populations around the world. The growth of cultured bivalves depends on food (e.g. phytoplankton) availability and quality as well as thermal conditions, which play an important role in physiological processes (Cardoso et al. 2007, Dutertre et al. 2009, Hollarsmith et al. 2020), with consequences for production and profitability. In addition, variations in density associated with farming practices or natural processes such as mortality may have feedback effects on production linked to local carrying capacity. Local food depletion induced by high-density bivalve filter-feeding activities can reduce the phytoplankton concentration in the water column, thereby affecting individual growth of the cultured bivalves (Ferreira et al. 2007, Duarte et al. 2008, Mazón-Suástegui et al. 2008, Cugier et al. 2010, Guyondet et al. 2010, Guyondet et al. 2013, Filgueira et al. 2014b, Bacher et

Market-size oysters, farming area and simulated oyster growth performance (blue: low, red: high) in Bourgneuf Bay.

Photos/Graphic: P. Cugier, Ifremer; SMIDAP, IGN 2022

© The authors 2022. Open Access under Creative Commons by Attribution Licence. Use, distribution and reproduction are unrestricted. Authors and original publication must be credited.

<footnote n="1">¹Ifremer, DYNECO, 29280 Plouzané, France</footnote>

<footnote n="2">²University of Brest, CNRS, IRD, Ifremer, LEMAR, 29280 Plouzané, France</footnote>

<footnote>*Corresponding author: philippe.cugier@ifremer.fr</footnote>
al. 2019). During episodes of high mortality, oyster farmers face drastic reductions in biomass, which can impact the remaining oysters' growth performance due to changes in trophic balances, ultimately affecting the quantity and quality of marketable oysters. Around the world in the last decades, oyster farmers have had to face episodes of massive oyster mortality linked to various pathogens (EFSA AHAW 2015). In many French farming sites, young oysters have been affected by a microvariant of the ostreid herpesvirus 1 (OsHV-1 µVar) and adult oysters by the bacteria *Vibrio aestuarianus*, which caused the death of more than 80% of the juveniles and 25% of the adults in the farms (Pernet et al. 2012, 2019, Gangnery et al. 2019).

Tools and management methods are needed to help bivalve farmers adapt to multiple sources of variation in production (Aubert et al. 2020). Different types of mathematical models have been developed to assess ecological carrying capacity (Weitzman & Filgueira 2020). Their complexity varies both in the accuracy of the processes, state variables taken into account and the way they represent space (Dabrowski et al. 2013). Several process-based models detail the interactions between bivalves and the ecosystem to assess the impact on both individual performance and ecosystem functioning or to project system changes. These models are widely applied in many mussel and oyster farming sites (Ferreira et al. 2008, Grangeré et al. 2009, 2010, Guyondet et al. 2010, Nunes et al. 2011, Dabrowski et al. 2013). Regarding ecophysiological models, many studies have applied the dynamic energy budget (DEB) theory (Kooijman 2010) to simulate growth and bio-energetics of individual bivalves and explain the spatial and temporal variability in their physiological performances (Bourlès et al. 2009, Rosland et al. 2009, Ren et al. 2010, Alunno-Bruscia et al. 2011, Barillé et al. 2011, Thomas et al. 2011). Some models have addressed food depletion in the water column due to bivalve filtration activity (Guyondet et al. 2005, 2013, Duarte et al. 2008, Cugier et al. 2010). Other models have dealt with ecological or shellfish carrying capacities (Duarte et al. 2003, Guyondet et al. 2010, 2015, Filgueira et al. 2014a). Several studies have assessed the impact of global changes, including those related to climate, on the growth and biogeography of cultivated bivalves (Guyondet et al. 2015, Filgueira et al. 2016, Thomas et al. 2016, 2018, Thomas & Bacher 2018, Steeves et al. 2018). Modelling tools have also been developed to assess disease propagation amongst individuals and between farms and

to potentially help the producers facing these phenomena (Thrush et al. 2017, Gangnery et al. 2019, Lupo et al. 2019, 2020). Due to their complexity, few models allow the combination of a large number of changing scenarios to help with decision making. This explains why models aimed to support production and management usually consider a limited number of components or have low spatial accuracy (Ferreira et al. 2007, Filgueira et al. 2014b). There is, therefore, a challenge in designing a simple modelling tool to predict realistic responses at spatial and temporal scales relevant for stakeholders.

On the French Atlantic coast, Bourgneuf Bay is a region of extensive oyster *Crassostrea gigas* aquaculture. This site exhibits strong spatio-temporal variability in food, turbidity and temperature that affects oyster growth performance (Palmer et al. 2020). One of the specificities of this sector is the strong contribution of microphytobenthos to primary production, mainly in the north region of the bay (Méléder et al. 2003, Barillé et al. 2011, Kazemipour et al. 2012). In this region, oyster farming has been facing mortality events for many years, similar to many other sites in France (Fleury et al. 2020a, Le Bihan et al. 2020). Producers often have to adapt quickly to these events. They mainly rely on their empirical knowledge built on past experiences to anticipate variations in production. The complexity of the interactions, spatial structuring and inter-annual variability, however, make it difficult to anticipate future variations in production. To support producers in the spatial planning of aquaculture activities which enable sustained production while maintaining control of the impacts on ecosystem production, effective management tools are needed.

This paper aims to assess the sensitivity of oyster production in Bourgneuf Bay to mortality, stocking density and variations in environmental conditions. Our work is based on the coupling of a 3D model of the bay ecosystem that simulates environmental conditions and an individual oyster growth model based on DEB theory. After validation of the model, the definition of a set of scenarios crossing various environmental, farming and mortality conditions yielded a total of 48 simulations that allowed evaluation, quantification and ranking of their impacts on oyster growth performance and production. We also aimed to build a simple tool based on this complex modelling approach to test management options. To this end, simulations were combined to produce a mathematical function that adapts rearing strategies to environmental variability and oyster mortality.
## 2. MATERIALS AND METHODS

### 2.1. Study area

Bourgneuf Bay is a semi-enclosed coastal area located along the Atlantic coast south of the Loire estuary (Fig. 1). It is a macrotidal area with a semi-diurnal tidal regime. The bay is largely open to the ocean in the northwest and by a narrow channel in the south. A large intertidal zone of 100 km² covers the eastern shore of the bay and is a region of extensive oyster *Crassostrea gigas* aquaculture, which occupies an area of ca. 1000 ha (Fig. 1). Standing stock is about 45 000 t, with an annual production of 10 000 t (Haure et al. 2003); we use these values as a baseline scenario for our simulations.

Based on a hydrodynamic model, Lazure (1992) established the main hydrodynamic characteristics of the bay. The strongest tidal currents occur in the open channel in the northern part of the bay (~1 m s⁻¹) and in the narrower, southern channel (~1 m s⁻¹). Inside the bay, velocities decrease to a maximum of 0.4–0.5 m s⁻¹. Residual circulation follows a general pattern, with a main transit from north to south and the presence of an eddy zone in the western part of the bay along the coast of Noirmoutier Island. In the eastern part of the bay, along the intertidal areas, residual currents become extremely weak. Nevertheless, residual circulation in the bay is highly variable, depending on wind conditions. The overall turnover time of water in the bay is about 2 mo. The water column is globally well-mixed thanks to the tidal currents and there is little to no stratification.

### 2.2. Modelling strategy

The modelling strategy is illustrated in the workflow diagram in Fig. 2. It is based on a 3D ecosystem model coupled with an individual oyster growth model based on the DEB theory. After a validation step, the model was applied, with a series of scenarios aiming to quantify the impact of stocking density, mortality and inter-annual environmental variability on the annual growth of oysters. Results were then used to establish an empir-

<figure><img src="image_6.png" /><figcaption>Fig. 1. Bourgneuf Bay, France, showing the distribution of cultivated oysters. The model was calibrated at Coupelasse (a) and Graisselous (b) stations; Coupelasse was also used for oyster growth validation</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 2. Workflow of the modelling steps, with references to sections and figures described in this paper</figcaption></figure>
ical response function that allows generalisation of the predictions. This empirical function allowed us to explore different management measures to achieve production objectives in response to inter-annual environmental variability and mortality. Each step and tool are detailed in the following sections.

## 2.3. Ecosystem model

Hydrodynamic characteristics were computed with the MARS 3D hydrodynamic model (Lazure & Dumas 2008). It uses a regular rectangular horizontal grid cell of 800 × 800 m (Fig. S1 in Supplement 1 at www.int-res.com/articles/suppl/q014p053_supp.pdf; for all supplements) with a vertical resolution of 10 sigma layers to compute the current velocity, temperature and salinity fields of Bourgneuf Bay. The Nutrient-Phytoplankton-Zooplankton-Detritus (NPZD) model, ECO-MARS3D, was coupled to hydrodynamics to simulate primary production. The main processes and equations of this module are described in Vahoutte-Brunier et al. (2008) and were revisited and recalibrated recently by Ménesguen et al. (2019) for the English Channel and Bay of Biscay. This model resolves the nutrient cycles of nitrogen (N), phosphorus and silica as well as the N content of 3 phytoplanktonic compartments: diatoms, dinoflagellates and nanoflagellates. Total chlorophyll a (chl a) is derived from planktonic nitrogenous state variables using an empirical chl a:N ratio, computed following a Smith-like formula (Smith 1936) to account for local light extinction.

In Bourgneuf Bay, microphytobenthic biomass is particularly important, especially in the northeastern part of the bay and has been studied and mapped for many years (Méléder et al. 2003, 2005, Kazemipour et al. 2012). Microphytobenthos play a crucial functional role in the bay, as they contribute to oyster diet (Cognie et al. 2001, Decottignies et al. 2007). Since the ECO-MARS3D model does not simulate microphytobenthos production, a flux of diatoms ($F_{mpb}$) was introduced at the bottom of the water column and spatially distributed on the basis of the known spatial distribution of microphytobenthos (Méléder et al. 2003, 2005, Kazemipour et al. 2012). $F_{mpb}$ is expressed in µmol N l⁻¹ d⁻¹. $F_{mpb}$ depends on season and is included in the model equation as follows:

$$ \begin{array}{l} \text{if } U_b > U_{bcrit} \text{ then } \frac{dC_{\text{diat}}}{dt} = F_{\text{mpb}} \\ \text{otherwise } \frac{dC_{\text{diat}}}{dt} = 0 \end{array} \qquad (1) $$

where $C_{\text{diat}}$ is the water concentration of diatoms (µmol l⁻¹ of N), $t$ is time (d), $U_b$ is the bottom current velocity and $U_{bcrit}$ is the critical velocity above which sediment resuspension occurs and the flux is activated in order to mimic resuspension of benthic diatoms. $F_{mpb}$ and $U_{bcrit}$ were calibrated to reproduce the realistic chl a levels observed at Coupelasse and Graisselous, located in the north and south of the bay, respectively (Figs. 1 & S2 in Supplement 2).

The northern part of the bay (Coupelasse sector) is characterized by high turbidity, with very muddy sediment in the intertidal zone. Further south, in the Graisselous sector, turbidity is lower and associated with a coarser, sandy-muddy sediment (Dutertre et al. 2009). Similarly, chl a and particulate organic matter concentrations are also higher in the north than in the south (Dutertre et al. 2009). Studies have also shown that the proportion of chl a related to resuspended microphytobenthos is important throughout the bay, higher in the north and in spring tides when currents are stronger (Haure & Baud 1995, Barille-Boyer et al. 1997). These observations guided the calibration of the $U_{bcrit}$ and $F_{mpb}$ parameters in an attempt to reproduce the chl a levels observed at Coupelasse and Graisselous as accurately as possible. This resulted in a lower value of $U_{bcrit}$ in the north, reflecting muddier and more easily mobilized sediment than in the south as well as a higher value of $F_{mpb}$ in the north than in the south, in accordance with the observations of chl a concentrations in the water column. The introduction of a difference in fluxes between summer and winter allowed us to account for the natural cycle of microphytobenthic production.

In the northern part of the bay, $U_{bcrit} = 0.05 \text{ m s}^{-1}$; $F_{mpb} = 100 \text{ µmol N l}^{-1} \text{ d}^{-1}$ from April to September, $F_{mpb} = 20 \text{ µmol N l}^{-1} \text{ d}^{-1}$ otherwise. In the southern part of the bay, $U_{bcrit} = 0.1 \text{ m s}^{-1}$; $F_{mpb} = 10 \text{ µmol N l}^{-1} \text{ d}^{-1}$ from April to September, $F_{mpb} = 2 \text{ µmol N l}^{-1} \text{ d}^{-1}$ otherwise.

## 2.4. Oyster growth model

Individual oyster growth was simulated using a DEB model (Kooijman 2010). The DEB model allows us to simulate, in a varying environment, the flow of energy through an organism; from its uptake to its use for structural growth, development, reproduction and maintenance. It offers a generic framework that covers the full life cycle, from embryo to juvenile to adult. The DEB model provides a quantitative prediction of mass and energy balances and was used here to simulate the physical length of oys-
ters. Applied to adult oysters, the equations and parameters are fully described in Thomas & Bacher (2018). Values of parameters specific to the present study are given in Table 1.

In the DEB model, ingestion rate was simulated by a Holling Type II functional response ($f \in [0,1]$) considering the particulate inorganic matter (PIM) concentration (Kooijman 2006) and modulated by a temperature correction function ($T_C$; see Eq. 3):

$$f = \frac{X}{X + K(Y)} \quad \text{with} \quad K(Y) = X_k \left(1 + \frac{Y}{X_{ky}}\right) \quad (2)$$

where $X$ is the food concentration (mg C ml$^{-1}$), $X_k$ the food half-saturation coefficient, $Y$ is PIM concentration (mg l$^{-1}$) and $X_{ky}$ is the PIM half-saturation coefficient. In previous applications of the DEB-oyster model, Grangeré et al. (2009, 2010) showed that oyster dry weight was better simulated using phytoplankton C concentration instead of chl a concentration as a quantifier for food. The chl a output of the model was converted to C through the empirical formulation of the chl a:C ratio given by Cloern et al. (1995) and was dependent on temperature, nutrient concentrations and irradiance. Since DEB model variables are expressed in energy (J), a conversion factor between energy and phytoplankton C is required. As done by Grangeré et al. (2009, 2010), a fixed ratio of 47.76 J mg$^{-1}$ C, calculated from Platt & Irwin (1973) was also used.

$T_C$ is formulated as follows:

$$T_C = \exp\left(\frac{T_A}{T_1} - \frac{T_A}{T}\right) \left[1 + \exp\left(\frac{T_{AL}}{T_1} - \frac{T_{AL}}{T_L}\right) + \exp\left(\frac{T_{AH}}{T_H} - \frac{T_{AH}}{T_1}\right)\right] \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]^{-1} \quad (3)$$

<table><caption>Table 1. Values of the dynamic energy budget (DEB) parameters used in the model to simulate individual oyster growth. PIM: particulate inorganic matter</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Food half saturation coefficient</td><td>X<sub>k</sub></td><td>310</td><td>µg C l<sup>-1</sup></td></tr><tr><td>PIM half saturation coefficient</td><td>X<sub>ky</sub></td><td>10</td><td>mg l<sup>-1</sup></td></tr><tr><td>Maximum surface specific ingestion rate</td><td>IR<sub>max</sub></td><td>1027</td><td>J cm<sup>-2</sup> d<sup>-1</sup></td></tr><tr><td>Arrhenius temperature</td><td>T<sub>A</sub></td><td>5800</td><td>K</td></tr><tr><td>Arrhenius temperature for lower boundary</td><td>T<sub>AL</sub></td><td>75000</td><td>K</td></tr><tr><td>Arrhenius temperature for upper boundary</td><td>T<sub>AH</sub></td><td>30000</td><td>K</td></tr><tr><td>Lower boundary tolerance range</td><td>T<sub>L</sub></td><td>281</td><td>K</td></tr><tr><td>Upper boundary tolerance range</td><td>T<sub>H</sub></td><td>300</td><td>K</td></tr></tbody></table>

where $T_A$ is the Arrhenius temperature within the tolerance range of the oysters (i.e. without physiological damages); $T_1$ is the temperature given $T_C = 1$; $T_L$ and $T_H$ the low and high boundaries of the tolerance range; and $T_{AL}$ and $T_{AH}$, the Arrhenius temperatures beyond the lower and higher boundaries.

The NPZD and oyster models are fully coupled, as the latter removes phytoplankton from the water column by the filtration process. The individual filtration rate (FR; J d$^{-1}$) is calculated as a function of the surface area of the organism (from structural volume, $V^{2/3}$), the area maximum specific ingestion rate (IRmax), temperature and food availability (Grangeré et al. 2009):

$$FR = IR_{\max} \frac{X}{(K[Y])} V^{2/3} T_c \quad (4)$$

In each grid cell of the model domain, an averaged FR for the entire population was estimated by multiplying FR (corresponding to an average individual) by the total number of individuals and then converting from J to chl a with the previously defined coefficients.

## 2.5. Forcing variables and initial conditions

The French Naval Hydrographic and Oceanographic Service (SHOM) provided the bathymetry of Bourgneuf Bay. The databases of the French water agencies were queried to obtain the main river inputs (average daily flows and nutrient concentrations based on monthly or bimonthly measurements). Outputs from the ARPEGE meteorological model of Météo-France (French weather forecasting service), with a spatial resolution of 0.5° and a temporal resolution of 6 h, provided meteorological data (air temperature, atmospheric pressure, air humidity, cloud cover, wind speed and direction) for forcing at the ocean-atmosphere interface. SeaWiFS, MODIS and MERIS satellite images were used to obtain daily concentrations of PIM. These images were merged as described in Saulquin et al. (2011) using a regional algorithm specifically designed for the coastal waters of the Bay of Biscay (Gohin 2011). They were interpolated on the computational grid and then used as inputs to the ECO-MARS3D model to calculate light extinction in the water acting on phytoplanktonic production, but also as a factor influencing the FR of oysters (see Eqs. 2 and 4).

The actual distribution of oyster densities was obtained from a standing stock assessment made by IFREMER in 2002 in Bourgneuf Bay (Haure et
al. 2003) and spatially distributed in grid cells with respect to farm locations and used in the model validation step. Finally, theoretical stocking density distributions were also used as described in Section 2.6.

A spin-up of 1 yr was run with the NPZD model in order to obtain stable environmental conditions before coupling with the oyster model. Details on model validation are given in Supplement 2.

## 2.6. Scenario building

The model was used to evaluate how the individual growth and local biomass of *C. gigas* are affected by variations in density and by interannual variation of hydroclimatic conditions described in the following sections.

### 2.6.1. Density scenarios

In Bourgneuf Bay, oysters are grown in bags placed on tables positioned about 1 m above the sea floor. The stock is divided into 2 age classes (cohorts) of 1 and 2 yr, respectively (Haure et al. 2003). For each density scenario, the standing stock was divided into 2 cohorts (C1 and C2) composed of individuals initially 6 and 18 mo old, respectively. A reference density (R) of 230 000 ind. ha⁻¹ was set for each cohort. This initial density allowed us to obtain a simulated average standing stock of approximately 40 000 t, which was close to the estimation of 46 000 t made in 2002 by Haure et al. (2003). We also considered 3 conditions covering a wide range of standing stock densities: (1) low density = R × 0.5 (115 000 ind. ha⁻¹); (2) high density = R × 2 (460 000 ind. ha⁻¹); and (3) very high density = R × 2.6 (600 000 ind. ha⁻¹). For all scenarios, the oyster stocks were homogeneously distributed across the farming area in the bay, which represented 52 model grid cells.

### 2.6.2. Mortality scenarios

We considered 2 mortality rate scenarios based on monthly observations made at Coupelasse station by the national shellfish observation network (RESCO) (Fleury et al. 2020b; Fig. S4 in Supplement 3): (1) a reference standard monthly mortality rate corresponding to a cumulative annual mortality of around 20 and 10% for cohorts C1 and C2, respectively; and (2) a high monthly mortality rate corresponding to a

cumulative annual mortality of 60 and 25% for cohorts C1 and C2, respectively.

We assumed that individual growth is not affected by the presence of a potential disease linked to mortality. Thus, these scenarios corresponded to a modulation of density that progressively decreases through the year with different kinetics depending on mortality rates.

### 2.6.3. Climatic conditions

We compared 2 years with contrasting hydroclimatic conditions. The year 2005 was considered a dry year, mainly characterized by a weak flow of the Loire River compared to the average flow estimated for the period 1950–2011 (see Fig. S5 in Supplement 3). Conversely, 2008 was considered a stormy and wet year, characterized by strong flows of Loire River in spring and autumn (Gohin et al. 2015).

### 2.6.4. Initial weight of C2 individuals

A first simulation was run for 1 yr with the initial weight of C1 and C2 based on averaged observations. We built scenarios of initial weight of C2 using the distribution of final weight of C1. The average weight is taken as the reference weight, with the 10th and 90th percentiles as low and high weight, respectively.

### 2.6.5. Simulation outputs

Simulations were run on an annual basis, from the 1st of January to the 31st of December. The growth performance of oysters (i.e. final weight) and total biomass in each cell were evaluated at the end of each simulated scenario.

Combinations of density, mortality, climatic conditions and initial weight scenarios were run for the year 2005 (considered the reference year; Table 2). Only few combinations were run for 2008, for comparison purposes. We defined the reference scenario as a combination of reference density, reference mortality, reference initial weight and year 2005 conditions.

Results of scenarios were interpreted by comparing the final weight obtained in each grid cell with the final weight of the reference scenario and by mapping the percentage of the difference, DIFᵢ, calculated as:

$$DIF_i = 100 (W_{H}^{S} - W_{H}^{R}) / W_{H}^{R} \quad (5)$$
<table><caption>Table 2. Simulated oyster *Crassostrea gigas* growth and production scenarios. Simulations were repeated for the reference and high scenarios of mortality, giving a total of 58 simulations. C1: cohort 1; C2: cohort 2; Lw: low initial weight; Rw: reference initial weight; Vw: very high initial weight; L: low density; H: high density; V: very high density</caption><thead><tr><th>Initial weight of C2</th><th>Initial density of C1</th><th>Initial density of C2</th><th>2005</th><th>2008</th></tr></thead><tbody><tr><td rowspan="3">Lw</td><td rowspan="2">L</td><td>L</td><td>X</td><td></td></tr><tr><td>V</td><td>X</td><td></td></tr><tr><td>V</td><td>L</td><td>X</td><td></td></tr><tr><td rowspan="12">Rw</td><td rowspan="3">L</td><td>L</td><td>X</td><td>X</td></tr><tr><td>R</td><td>X</td><td></td></tr><tr><td>H</td><td>X</td><td></td></tr><tr><td rowspan="3">V</td><td>L</td><td>X</td><td>X</td></tr><tr><td>R</td><td>X</td><td>X</td></tr><tr><td>H</td><td>X</td><td></td></tr><tr><td rowspan="3">H</td><td>L</td><td>X</td><td></td></tr><tr><td>R</td><td>X</td><td></td></tr><tr><td>H</td><td>X</td><td></td></tr><tr><td rowspan="3">V</td><td>L</td><td>X</td><td></td></tr><tr><td>R</td><td>X</td><td></td></tr><tr><td>H</td><td>X</td><td></td></tr><tr><td rowspan="6">Vw</td><td rowspan="2">L</td><td>L</td><td>X</td><td></td></tr><tr><td>V</td><td>X</td><td></td></tr><tr><td rowspan="2">V</td><td>L</td><td>X</td><td></td></tr><tr><td>V</td><td>X</td><td></td></tr></tbody></table>

where $W_{i1}^s$ is the final weight in cell $i$ in scenario s and $W_{i1}^R$ is the final weight in cell $i$ in the reference scenario. In order to compare scenarios, a mean sensitivity index (SI) was computed as:

$$ SI = \frac{1}{n} \sum_{i=1,n} \text{DIF}_i \qquad (6) $$

where $n$ is the total number of cells.

## 2.7. Response function deduced from simulated scenarios

### 2.7.1. Function building

Our coupled 3D biogeochemical and oyster growth models allowed us to map the responses in terms of individual growth and production for a range of environmental changes, but there were limits to computing time and thus the number of simulations and the

opportunity to examine new scenarios was limited. Therefore, we combined the outputs of all of our simulations and built a simple response function that predicts the final weights of the 2 cohorts in each cell (see equations in Supplement 4). The function uses initial numbers of oysters, initial oyster weights and the mortality rates for the 2 cohorts as input data. We obtained as many functions as there were cells and years (2005 or 2008). Function outputs were the predicted number and weight of oysters in C1 and C2 at the end of the year. The annual production was defined as the product of the number and the weight of oysters in C2 at the end of the year

### 2.7.2. Management measures

We considered management measures based on the seeding of C1 individuals. Our reference scenario corresponded to low mortality rates and an environment similar to that of 2005. We applied the response function to a series of seedings with the same mortality rates and environment. We summed the values over all cells to assess annual production over the cultivated area. Production was then plotted as a function of seeding to show the relationship between seeding and production. We also plotted average final weight of C2 oysters as a function of seeding.

We then defined a new scenario with a high mortality rate for the 2 cohorts and repeated the same calculations for a series of seedings. We compared this new scenario to the reference scenario by plotting 2 curves showing the relationship between seeding and production and did the same for the 2 curves linking seeding and the final weight of C2 oysters. We assumed that farmers would increase seedings to reach the same level of production in the scenario with high mortality as in the reference scenario. Using the curve plots derived from the series of calculations, we assessed how much seeding should be increased to achieve the same production or final weight as the reference scenario.

We also compared the same reference scenario with a scenario based on the same mortality rates but a different environment. In this case, we repeated the calculations using the parameters of the response curve based on 2008 scenarios and again plotted the annual production and the final individual weight versus seeding.

Finally, we considered a series of scenarios based on a combination of mortality rates for the 2 cohorts. Using an optimization algorithm, in each case, we calculated the optimal seeding which would allow
the farmers to achieve the same level of production as the reference production. Our objective function was the absolute difference between production in the reference scenario and the scenario with new mortalities and depends only on seeding. All estimations were computed with the ‘optimize’ function in R (R Core Team 2020). Data collected through the monitoring of oyster mortality in Bourgneuf Bay showed that C1 mortality varied from <0.1 to 0.8 and that C2 mortality varied from ca. 0.05 to 0.25. Using these ranges of mortality, we displayed the results as a surface curve showing optimum seeding for any combination of mortality rates.

## 3. RESULTS

### 3.1. Spatial patterns of growth

Under reference conditions, simulated growth performance appeared strongly heterogeneous throughout the rearing area and exhibited similar patterns for the 2 cohorts (Fig. 3A). Higher performances were found in the northern and western parts of the bay, with final individual total weights reaching 40 and 90 g for C1 and C2, respectively. Lower performances were seen in the southern part of the bay, with final total weights < 15 g for C1 and 40 g for C2.

This spatial heterogeneity of growth performance was associated with strong heterogeneity in environmental conditions in the bay (Fig. 3B,C,D). These conditions modulate the energy acquisition capacities of the oysters, here illustrated by the average functional response in the rearing areas (Fig. 3E).

### 3.2. Inter-annual variability

Simulations showed strong differences in final weights between the reference year (2005) and the stormy and wetter year (2008), with an average (±SD) reduction of -29 ± 6 and -18 ± 5% for C1 and C2, respectively (Fig. 4). These changes were not homogeneous across the bay. The western and central parts of the bay were more sensitive, with a decrease in growth performance of -39 and -31% for C1 and C2, respectively.

The simulated interannual variations in growth potential and associated spatial heterogeneity were mostly sustained by differences in feeding rates. The mean annual functional response (f) anomaly is always negative when mean annual T~C~ exhibits both positive and negative values, depending on locations

(Fig. 5A). However, the average anomalies remain at low values. Significant differences between the 2 simulated years are found when the cumulative values of f and T~C~ are represented (Fig. 5B). The comparison of cumulated T~C~ does not differ between 2005 and 2008, whereas a clear difference is observed for cumulated f. For 2008, the cumulated f remains close to zero until the end of April, reflecting very low feeding rates, whereas for 2005 the cumulated f increases rapidly from March. This delay in food acquisition is never compensated for, which results in a net deficit for 2008, comparable to the observations made on the final biomasses. A strong difference in the timing of the spring bloom was observed between the 2 years, with a higher and earlier peak in 2005 compared to 2008, which explains the direct effect on oyster feeding potential. This difference in phytoplankton growth can be explained by light limitation due to suspended matter load. The beginning of 2008 was characterized by storm events, high river discharge of the Loire River, located north of Bourgneuf Bay (see Section 2.6.3), bringing higher turbidity to the inner bay. The higher turbidity may also partly explain the differences in growth, as it limits the food intake of the oyster. Thus, spring bloom is a critical period for oysters and supports a significant portion of the annual growth. This is illustrated in Fig. S7 in Supplement 5, which shows sea-surface suspended matter measured by satellite and averaged from January to March as well as the time series of observed suspended matter and simulated chl a averaged over the bay for 2005 and 2008.

### 3.3. Response of final weight to initial density level

Increased initial densities led to a general decrease in final weight for the 2 cohorts, with higher reductions in lower food concentration areas located in the south-eastern part of the bay (Fig. 6). C1 was more affected than C2, with a final weight reduction between -15 and -35% for C1 and -5 and -20% for C2.

### 3.4. Response of final weight to mortality levels

Increased mortality reduced the density of the oysters, which had a positive effect on individual growth. An increase of 4–8% and 0–3% of the final weight was seen for C1 and C2, respectively (Fig. 7). This effect was greatest in the south-eastern part of the bay, where growth performance is naturally lower.
<figure><img src="image_8.png" /><figcaption>Fig. 3. Oyster growth patterns in Bourgneuf Bay. (A) Oyster final weight simulated in each model cell for the 2 cohorts: C1 and C2. Results are shown for 2005 and the reference scenarios of density, initial weight of C2 and mortality. (C–D) Annual averages of the 2 variables simulated by the model: (C) chl *a*, (D) temperature. (E) Annual average of suspended particulate inorganic matter (SPIM) observed by satellite and used as forcing in the model. Black dots: cells with oyster rearing activity. (F) Yearly mean functional response (f) modulated by the temperature correction factor ($T_C$), applied in the DEB model</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Fig. 4. Inter-annual variability of oyster final weight expressed as a percentage of the difference between scenarios in 2008 and 2005. Results are given for the 2 cohorts (C1 and C2) for reference initial densities and reference mortality of C1 and C2</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 5. Inter-annual variation in forcing variables. (A,B) Mean annual anomaly computed after the (A) daily-scaled functional response (f) and (B) temperature correction factor (T<sub>C</sub>). (C,D) Time series of cumulated (C) f and (D) T<sub>C</sub> for the scenarios in 2005 and 2008; Thin lines: single cells; thick lines: mean over all cells</figcaption></figure>
### 3.5. Response of final weight to initial weight of C2

Fig. 8 illustrates the effect of the initial weight of individuals in C2 (see Section 2.6.4) on the annual growth rate of C1 and C2 (expressed in g yr⁻¹) depending on whether the initial density of C1 and C2 were low or very high. In this example, the mortality was chosen as the reference mortality. As the initial weight of individuals in C2 increases, the growth rate of individuals in C1 decreases with a greater effect for the high initial density of C2. As the biomass of C2 increases (due to density or individual weight), trophic competition increases, leading to less available resources for C1 individuals.

For C2, annual growth tends to increase slightly but quickly reaches a maximum when the initial weight of C2 increases. On the other hand, there was almost no gain in growth, or even a decrease, when the initial density of the C2 cohort was very high, as the reared biomass has probably exceeded the trophic capacity.

### 3.6. Sensitivity analysis

SI was used to compare scenarios (Fig. 9) and ranged from -30 to +30%. The highest SI resulted from the interannual comparison between the 'wet' year (2008) and the 'dry' year (2005). The SI reached -30% for C1 and -20% for C2. An increase in C2 density had the second-largest impact and led to a decrease of 12 and 18% in the final weight of C2 and

C1, respectively. An increase in the initial mass of C2 individuals is a special case because this intrinsically has a strong impact on the final weight of C2 individuals (>20%), especially as the annual growth rate increases with the initial weight (see the previous section). The smallest variations in SI (<10%) were related to mortality and density variations.

### 3.7. Scenarios of management measures

We used the response function to compute and plot annual production and final individual weight for a series of seedings. Calculations clearly showed the inverse relationship between the final individual weight and seedings and the non-linear increase of production (Fig. 10, blue lines for reference scenario). We then compared these relationships for the reference scenario and a scenario with high mortality (Fig. 10, red lines). Increasing mortality decreases oyster density and improves individual growth, which explains the upwards shift of the curve relating seeding with final weight. However, production with high mortality always remains below production with reference mortality, which means that the gain in individual weight is not sufficient to compensate for losses due to mortality. The relationship between seeding and production allows us to infer that seeding equal to ca. 5.6 × 10⁵ ind. ha⁻¹ would allow us to retrieve the production of the reference scenario corresponding to a seeding of ca. 2.3 × 10⁵ (Fig. 10, vertical lines).

<figure><img src="image_11.png" /><figcaption>Fig. 6. Effect of initial density scenarios on oyster growth. Oyster final weight variation was computed according to the high and reference density scenarios for cohorts C1 and C2. Results are shown for year 2005, reference initial weight of C2 and reference mortality</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 7. Effect of mortality scenarios on oyster growth. Oyster final weight variation is expressed as percentage variation between the reference and high mortality scenarios for cohorts C1 and C2. Results are given for year 2005, reference initial density of C1 and C2 and reference initial weight of C2</figcaption></figure>

Similar patterns emerged when we carried out the same analysis using different environmental conditions, i.e. environment corresponding to 2008, without any change in mortality rate ind. ha⁻¹ (Fig. 10, green lines). Increasing density to ca. 3.9 × 10⁵ ind. ha⁻¹ would allow us to achieve the production of the reference scenario, but the final weight would remain much lower, i.e. ca. 25 g instead of ca. 40 g as in the reference scenario.

The positive relationship between seeding and production makes it possible to estimate an exact seeding rate, allowing us to achieve the same production as the reference scenario when mortality rates of C1 and C2 vary. Combining sets of mortality rates, we estimated that seeding would vary from 2 × 10⁵ to 9 × 10⁵ ind. ha⁻¹ (Fig. 11). Farmers could use this rule, based on a new seeding, to precisely compensate for mortality loss. Similar calculations can

<figure><img src="image_13.png" /><figcaption>Fig. 8. Effect of initial weight of oyster cohort C2 on the individual growth rate of cohort C1 and C2. Two initial density scenarios are compared: low density of C1 and C2 and low density of C1 and very high density of C2</figcaption></figure>
<figure><img src="image_14.png" /></figure>

Fig. 9. Sensitivity of oyster growth to each factor tested (mean [± SD] spatial variability of final weight expressed as the percentage of difference compared to the reference). Results of the 5 factors tested are shown for the 2 cohorts (C1 and C2). We distinguished the effect of C2 initial weight on C2 final weight (hatched bar) which reflects the effect of a variable on itself

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

Fig. 10. (A) Production (t ha⁻¹) and (B) final individual weight (g) of cohort C2 computed for a series of seedings (ind. ha⁻¹) with an increase in mortality or change in environmental conditions. Blue line: reference mortality and year 2005; red line: increased mortality and year 2005; green line: reference mortality and year 2008. Vertical blue line: seeding value in the reference scenario (ca. 2.3 × 10⁵ ind. ha⁻¹); vertical red and green lines: seedings required to achieve the same level of production as the reference scenario. Dotted lines: changes in seedings needed to maintain production

easily be implemented in other scenarios where the growth rate is modified by environmental factors (e.g. 2008). Not shown here but easy to implement with our response function, another management measure could consist of increasing the initial weight of C1 individuals at seeding; this would modify the time needed to reach marketable size.

# 4. DISCUSSION

## 4.1. Drivers of oyster growth

Our study allows a better understanding of the spatial variability of oyster growth in Bourgneuf Bay. The best growth performances corresponded to areas where accessibility to food is optimal. Growth performance depends on the absolute availability of chlorophyll (Ren & Ross 2001, Ren et al. 2010) but also on temperature and the concentration of suspended matter (Bourlès et al. 2009, Barillé et al. 2011), which influence physiology and physical assimilation processes. Since spatial gradients of temperature are relatively weak, the analysis of the simulated functional responses clearly demonstrates that growth patterns result from the combined gradients of chl *a* and PIM. Simulations predicted that the best growth performances are in the northern and western parts of the bay. In the north, despite high turbidity linked to strong resuspension, oyster growth is achieved due to high concentrations of chl *a*, and it is likely that the growth performance would be higher if turbidity was lower. The high levels of PIM in the north certainly limit the capacity of bivalves to assimilate the available food and thus restrict growth in that area. Conversely, in the western part of the bay, lower concentrations of chl *a* but very low levels of turbidity result in good growth conditions comparable to those in the north. Spatial gradients of chl *a* are characterised by high concentrations in the north part of the bay linked to important microphytobenthic production on the wide foreshore of this area (Méléder et al. 2003, 2005, Kazemipour et al. 2012), which has been highlighted as a major component of the oyster diet (Cognie et al. 2001, Decottignies et al. 2007). Incorporating the microphytobenthos in an ecosystem model, for the first time even empirically, enabled a realistic simulation of the chlorophyll gradient (see validation in Supplement 2) with a direct impact on the simulated growth of oysters.

Food depletion is a feedback mechanism that can also affect food concentration and individual growth. At high rearing densities, filtration pressure
<figure><img src="image_16.png" /><figcaption>Fig. 11. Value of seeding cohort C1 (×10⁵ ind. ha⁻¹) that would allow us to achieve the same level of production as the reference scenario when mortality of C1 (x-axis) or C2 (y-axis) varies. Dots: observed mortalities reported in the RESCO monitoring network (Fleury et al. 2020b)</figcaption></figure>

decreases food concentration in the water column and, in extreme conditions, leads to trophic competition among individuals. For instance, Cugier et al. (2010) showed that cultivated filter feeders (oysters and mussels) can reduce the annual maximum chlorophyll locally by half in Mont-Saint-Michel Bay. They simulated scenarios with oyster and mussel rearing densities and showed that individual final weight could vary by several tens of percentage points depending on the scenario (Bacher et al. 2019). In another macrotidal bay, Grangeré et al. (2010) predicted that depletion due to filtration could reduce the final weight 2–4 times, depending on the area. In our study, simulations showed that the final weight of individuals in the 2 age classes decreased sharply, with a stronger effect on C1, when the density of C2 was doubled. The most affected areas are those in the south, close to the coast, where oyster growth is naturally lower. The rearing density is clearly too high compared to the available resources, and thus the system exceeds its carrying capacity. The spatial distribution of stocks in a given area can therefore be an important key to manage the growth of bivalves. However, we only considered that oyster stocks and changes in rearing densities were homogeneously distributed. Even though hydrodynamic studies have shown weak residual circulations on the foreshores (Lazure 1992), more complex scenarios could therefore be useful to assess how local changes in rearing density affect growth in other locations.

Our study clearly shows that interannual variability of hydroclimatic conditions is a major source of oyster growth variability. In contrasting situations, these variations can be much more important than those related to cultivation practices. From one year to the

next, climate conditions can be very different and can have a significant impact on key environmental variables for oysters growth: phytoplankton, temperature and inorganic material load. In coastal areas, the temporal variability of phytoplankton is usually correlated with meteorological conditions which increase nutrient inputs from catchment areas during rainy episodes and/or increase local turbidity during storm events. Although the maximum chlorophyll values reached in 2005 and 2008 in Bourgneuf Bay were very similar, the bloom started much later in 2008 than in 2005. Higher turbidity delayed the start of the spring growth of oysters, which was not compensated for during the rest of the year. In the English Channel, Grangeré et al. (2009) also showed that interannual fluctuations of oyster weight were consistent with variations in food availability (phytoplankton) and therefore with fluctuations in food assimilation. They also simulated a strong difference in growth between years, particularly linked to time shifts in spring blooms. Moreover, they showed that oyster growth also depends on the physiological state of the oysters at the beginning of spring (i.e. the status of their reserve) and therefore on the environmental conditions they experienced the previous fall and/or winter. More generally, the intensity of the oysters' metabolic activity during this period enhances or lowers their vulnerability during the rest of the year (Thomas et al. 2018). Thus, a low level of reserves in spring combined with a late phytoplankton bloom can further reduce oyster growth. We propose that carrying out simulations on 2 consecutive years with contrasting hydroclimatic conditions would provide a perspective to better evaluate the impact of the physiological status of oysters in spring on their growth later in the year.
## 4.2. Model use for management

Even if process-based models are useful for assess-
ing the main drivers of trophic capacity (Guyondet et
al. 2010, Nunes et al. 2011, Dabrowski et al. 2013),
extensive computing time and the need to use main-
frame computers limits the number of scenarios and
makes the application of such models to management
problematic. Various strategies can be implemented
to address this issue: (1) reduce the spatial dimensions
(e.g. hypotheses on horizontal and/or vertical homo-
geneity); (2) reduce the scale of the modelled out-
puts, e.g. bivalve individual scale vs. overall produc-
tion of the area (Ferreira et al. 2007); (3) simplify
model implementation by removing some feedback
processes (Tissot et al. 2012). In our study, we chose a
hybrid approach, starting with a complex 3D model
and accounting for the interactions between environ-
mental conditions and oyster growth. Based on this
validated model, we built a database of model pre-
dictions combining a range of input data, e.g. rearing
densities, mortality and environmental conditions.
Finally, these results allowed us to generate a simpli-
fied model in the form of response curves that could
be used to test different rearing strategies. The re-
sponse function we developed can help producers to
implement the best strategy to achieve given objec-
tives and/or face events affecting production. Appli-
cation of the response function simulating a loss of
production due to poor environmental conditions or
high mortality provided an estimate of the number of
C1 individuals required to be seeded in order to
return to initial (or reference) production. However,
other strategies can be envisaged, and the response
function, which is easy to use, can allow us to explore
them in association with farmers in the frame of com-
panion modelling approaches. In the case where
poor environmental conditions would result in a loss
of production, increasing the rearing density of C1 by
40% would allow a return to the reference produc-
tion level—but with a drastic reduction in the aver-
age individual final weight, from 44 to 25 g after 1 yr
of growth. Selling smaller oysters with a lower com-
mercial value may not be without financial conse-
quences for producers and the profitability of their
farms. Thus, in order to allow the oysters to reach a
minimum commercial weight, producers may have to
extend the rearing cycle, with a higher risk of haz-
ards linked to environmental or sanitary conditions
as well as consequences associated with the man-
agement of successive cohorts as a result of longer
immobilisation of stocks on the production sites.
Nunes et al. (2011) highlighted this type of effect and

showed that total mussel production increased con-
tinuously when the stocking density was gradually
increased, but that marketable size was reached
much later. In practice, producers would favour
lower densities and productions that allow them to
obtain larger mussels associated with a better market
price and a shorter growth cycle with lower eco-
nomic risk. The response function we developed can
thus be used to assess the conditions for achieving
more complex objectives accounting for the con-
straints of producers (e.g. minimum marketing size,
maximum length of rearing cycles).

Our hybrid approach addresses the need for adap-
tive management of aquaculture. Shellfish farmers
need to adapt their rearing strategy to technical,
administrative, economic and environmental con-
straints (Brigolin et al. 2017, Barillé et al. 2020). They
also face unpredicted risks such as episodes of
pathogen-related mortality which can lead to sharp
reductions in stocking densities, strongly affecting
their production and impacting their economic out-
come. Producers generally anticipate potential losses
by controlling spat collection, increasing the rearing
density of the young cohorts (Le Bihan et al. 2020),
but this practice is not without consequences for the
growth performance of the bivalves. Such a strategy
eventually results in a balance between the anticipa-
tion of mortality and the densities to be maintained in
rearing, which producers manage as well as possible
based on their empirical knowledge of the environ-
ment, the local and/or national dimension of the haz-
ard and the different forms of aquaculture practices
(Le Bihan et al. 2020). These choices, if made incor-
rectly, can limit the profitability of farms, i.e. the ben-
efits that producers can draw from their farms. These
benefits depend on the value of their production (the
market price of their product) but also on production
costs. Modelling approaches can help determine the
best strategy by integrating biological, epidemiolog-
ical and socio-economic components. Models allow
us to simulate how management measures that mod-
ify aquacultural biological performance (changes
in production areas, changes in stocking densities,
etc.) or external events (weather conditions, climate
change, mortality, etc.) could impact the economic
performance of farms or market prices (Merino et al.
2010, Kankainen et al. 2012). Our response function
approach for testing management measures on bio-
logical aquaculture performance is flexible and fully
capable of producing input for a bio-economic model
designed to assess the costs and benefits of manage-
ment measures and natural variabilities (mortality,
weather conditions).
*Acknowledgements*. This work received financial support from the National Research Agency (ANR) within the framework of the GIGASSAT (ANR-12-AGRO-0001) project.

LITERATURE CITED

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

Alunno-Bruscia M, Bourles Y, Maurer D, Robert S and others (2011) A single bio-energetics growth and reproduction model for the oyster *Crassostrea gigas* in six Atlantic ecosystems. *J Sea Res* 66:340–348

Aubert A, Aschenbroich A, Gaertner JC, Latchere O, Archambault P, Gaertner-Mazouni N (2020) Assessment of carrying capacity for bivalve mariculture in subtropical and tropical regions: the need for tailored management tools and guidelines. *Rev Aquacult* 12: 1721–1735

Bacher C, Gangnery A, Cugier P, Mongruel R, Strand Ø, Frangoudes K (2019) Spatial, ecological and social dimensions of assessments for bivalve farming management. In: Smaal AC, Gerreira JG, Grant J, Petersen JK, Strand Ø (eds) Goods and services of marine bivalves. Springer International Publishing, Cham, p 527–549

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

Barillé L, Lerouxel A, Dutertre M, Haure J, Barillé AL, Pouvreau S, Alunno-Bruscia M (2011) Growth of the Pacific oyster (*Crassostrea gigas*) in a high-turbidity environment: comparison of model simulations based on scope for growth and dynamic energy budgets. *J Sea Res* 66: 392–402

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

Barillé L, Le Bris A, Gouletquer P, Thomas Y and others (2020) Biological, socio-economic, and administrative opportunities and challenges to moving aquaculture offshore for small French oyster-farming companies. Aquaculture 521:735045

Barille-Boyer AL, Haure J, Baud JP (1997) L'ostréiculture en baie de Bourgneuf. Relation entre la croissance des huitres *Crassostrea gigas* et le milieu naturel: synthèse de 1986 a 1995. DRV/RA/RST/97-16. https://archimer.ifremer.fr/doc/00000/1633/

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

Bourlès Y, Alunno-Bruscia M, Pouvreau S, Tollu G and others (2009) Modelling growth and reproduction of the Pacific oyster *Crassostrea gigas*: advances in the oyster-DEB model through application to a coastal pond. *J Sea Res* 62:62–71

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

Brigolin D, Porporato EMD, Prioli G, Pastres R (2017) Making space for shellfish farming along the Adriatic coast. *ICES J Mar Sci* 74:1540–1551

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

Cardoso JFMF, Langlet D, Loff JF, Martins AR, Witte JI, Santos PT, van der Veer HW (2007) Spatial variability in growth and reproduction of the Pacific oyster *Crassostrea gigas* (Thunberg, 1793) along the west European coast. *J Sea Res* 57:303–315

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

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

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

Cognie B, Barillé L, Rincé Y (2001) Selective feeding of the oyster *Crassostrea gigas* fed on a natural microphytobenthos assemblage. *Estuaries* 24:126–134

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

Cugier P, Struski C, Blanchard M, Mazurié J and others (2010) Assessing the role of benthic filter feeders on phytoplankton production in a shellfish farming site: Mont Saint Michel Bay, France. *J Mar Syst* 82:21–34

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

Dabrowski T, Lyons K, Curé M, Berry A, Nolan G (2013) Numerical modelling of spatio-temporal variability of

growth of *Mytilus edulis* (L.) and influence of its cultivation on ecosystem functioning. *J Sea Res* 76:5–21

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

Decottignies P, Beninger PG, Rincé Y, Robins RJ, Riera P (2007) Exploitation of natural food sources by two sympatric, invasive suspension-feeders: *Crassostrea gigas* and *Crepidula fornicata*. *Mar Ecol Prog Ser* 334:179–192

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

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 Modell* 168:109–143

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

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

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

Dutertre M, Beninger PG, Barillé L, Papin M, Rosa P, Barillé AL, Haure J (2009) Temperature and seston quantity and quality effects on field reproduction of farmed oysters, *Crassostrea gigas*, in Bourgneuf Bay, France. *Aquat Living Resour* 22:319–329

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

EFSA AHAW (EFSA Panel on Animal Health and Welfare) (2015) Oyster mortality. EFSA J 13:4122

FAO (2018) The state of world fisheries and aquaculture 2018—meeting the sustainable development goals. FAO, Rome

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

Ferreira JG, Hawkins AJS, Bricker SB (2007) Management of productivity, environmental effects and profitability of shellfish aquaculture—the Farm Aquaculture Resource Management (FARM) model. *Aquaculture* 264:160–174

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

Ferreira JG, Hawkins AJS, Monteiro P, Moore H and others (2008) Integrated assessment of ecosystem-scale carrying capacity in shellfish growing areas. *Aquaculture* 275: 138–151

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

Filgueira R, Grant J, Strand Ø (2014a) Implementation of marine spatial planning in shellfish aquaculture management: modeling studies in a Norwegian fjord. *Ecol Appl* 24:832–843

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

Filgueira R, Guyondet T, Comeau LA, Grant J (2014b) A fully-spatial ecosystem-DEB model of oyster (*Crassostrea virginica*) carrying capacity in the Richibucto Estuary, Eastern Canada. *J Mar Syst* 136:42–54

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

Filgueira R, Guyondet T, Comeau LA, Tremblay R (2016) Bivalve aquaculture–environment interactions in the context of climate change. *Glob Change Biol* 22: 3901–3913

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

Fleury E, Barbier P, Petton B, Normand J and others (2020a) Latitudinal drivers of oyster mortality: deciphering host, pathogen and environmental risk factors. *Sci Rep* 10: 7264

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

Fleury E, Normand J, Lamoureux A, Bouget JF and others (2020b) RESCO REMORA database: national monitoring network of mortality and growth rates of the sentinel oyster *Crassostrea gigas*. *SEANOE*.

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

Gangnery A, Normand J, Duval C, Cugier P and others (2019) Connectivities with shellfish farms and channel rivers are associated with mortality risk in oysters. *Aquacult Environ Interact* 11:493–506

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

Gohin F (2011) Annual cycles of chlorophyll-a, non-algal suspended particulate matter, and turbidity observed from space and in situ in coastal waters. *Ocean Sci* 7:705–732

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

Gohin F, Bryère P, Griffiths JW (2015) The exceptional surface turbidity of the North-West European shelf seas during the stormy 2013–2014 winter: Consequences for the initiation of the phytoplankton blooms? *J Mar Syst* 148: 70–85

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

Grangeré K, Ménesguen A, Lefebvre S, Bacher C, Pouvreau S (2009) Modelling the influence of environmental fac-
tors 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

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

Grangeré K, Lefebvre S, Bacher C, Cugier P, Ménesguen A (2010) Modelling the spatial heterogeneity of ecological processes in an intertidal estuarine bay: dynamic interactions between bivalves and phytoplankton. Mar Ecol Prog Ser 415:141-158

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

Guyondet T, Koutitonsky VG, Roy S (2005) Effects of water renewal estimates on the oyster aquaculture potential of an inshore area. J Mar Syst 58:35-51

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

Guyondet T, Roy S, Koutitonsky VG, Grant J, Tita G (2010) Integrating multiple spatial scales in the carrying capacity assessment of a coastal ecosystem for bivalve aquaculture. J Sea Res 64:341-359

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

Guyondet T, Sonier R, Comeau LA (2013) Spatially explicit seston depletion index to optimize shellfish culture. Aquacult Environ Interact 4:175-186

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

Guyondet T, Comeau LA, Bacher C, Grant J, Rosland R, Sonier R, Filgueira R (2015) Climate change influences carrying capacity in a coastal embayment dedicated to shellfish aquaculture. Estuaries Coasts 38:1593-1618

Haure J, Baud JP (1995) Approche de la capacité trophique dans le bassin ostréicole (Baie de Bourgneuf). RIDRV-95-16/RA-BOUIN. https://archimer.ifremer.fr/doc/00000/6441/

Haure J, Martin JLY, Dupuy B, Nurry M and others (2003) Estimation des stocks d'huîtres en élevage dans la Baie de Bourgneufen 2002. Rapp IFREMER Novembre 2003 DRV/RA/LCPL/2003-4:20

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

Hollarsmith JA, Sadowski JS, Picard MMM, Cheng B, Farlin J, Russell A, Grosholz ED (2020) Effects of seasonal upwelling and runoff on water chemistry and growth and survival of native and commercial oysters. Limnol Oceanogr 65:224-235

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

Kankainen M, Setälä J, Berrill IK, Ruohonen K, Noble C, Schneider O (2012) How to measure the economic impacts of changes in growth, feed efficiency and survival in aquaculture. Aquac Econ Manag 16:341-364

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

Kazemipour F, Launeau P, Méléder V (2012) Microphytobenthos biomass mapping using the optical model of diatom biofilms: application to hyperspectral images of Bourgneuf Bay. Remote Sens Environ 127:1-13

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

Kooijman SALM (2006) Pseudo-faeces production in bivalves. J Sea Res 56:103-106

Kooijman SALM (2010) Dynamic energy budget theory for metabolic organisation, 3rd edn. Cambridge University Press, Cambridge

Lazure P (1992) Etude de l'hydrodynamique de la Baie de Bourgneuf. Ifremer report DEL-92.24, Ifremer's institutional repository. https://archimer.ifremer.fr

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

Lazure P, Dumas F (2008) An external–internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS). Adv Water Resour 31:233-250

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

Le Bihan V, Catalo M, Le Bihan J (2020) Reorganization of the value chain activities of oyster companies on the Atlantic coast following health crises in France (2006–2013). Mar Policy 117:103154

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

Lupo C, Travers MA, Tourbiez D, Barthélémy CF, Beaunée G, Ezanno P (2019) Modeling the transmission of *Vibrio aestuarianus* in Pacific oysters using experimental infection data. Front Vet Sci 6:142

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

Lupo C, Dutta BL, Petton S, Ezanno P and others (2020) Spatial epidemiological modelling of infection by *Vibrio aestuarianus* shows that connectivity and temperature control oyster mortality. Aquacult Environ Interact 12:511-527

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

Mazón-Suástegui JM, Ruíz-Ruíz KM, Parres-Haro A, Saucedo PE (2008) Combined effects of diet and stocking density on growth and biochemical composition of spat of the Cortez oyster *Crassostrea corteziensis* at the hatchery. Aquaculture 284:98-105

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

Méléder V, Launeau P, Barillé L, Rincé Y (2003) Cartographie des peuplements dus microphytobenthos par télé-détection spatiale visible-infrarouge dans un écosystème conchylicole. C R Biol 326:377-389

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

Méléder V, Barillé L, Rincé Y, Morançais M, Rosa P, Gaudin P (2005) Spatio-temporal changes in microphytobenthos structure analysed by pigment composition in a macrotidal flat (Bourgneuf Bay, France). Mar Ecol Prog Ser 297: 83-99

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

Ménesguen A, Dussauze M, Dumas F, Thouvenin B, Garnier V, Lecornu F, Répécaud M (2019) Ecological model of the Bay of Biscay and English Channel shelf for environmental status assessment. Part 1: nutrients, phytoplankton and oxygen. Ocean Model 133:56-78

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

Merino G, Barange M, Mullon C, Rodwell L (2010) Impacts of global environmental change and aquaculture expansion on marine ecosystems. Glob Environ Change 20: 586-596

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

Nunes JP, Ferreira JG, Bricker SB, O'Loan B and others (2011) Towards an ecosystem approach to aquaculture: assessment of sustainable shellfish cultivation at different scales of space, time and complexity. Aquaculture 315:369-383

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

Palmer SCJ, Gernez PM, Thomas Y, Simis S, Miller PI, Glize P, Barillé L (2020) Remote sensing-driven Pacific oyster (*Crassostrea gigas*) growth modeling to inform offshore aquaculture site selection. Front Mar Sci 6:802

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

Pernet F, Barret J, Le Gall P, Corporeau C and others (2012) Mass mortalities of Pacific oysters *Crassostrea gigas* reflect infectious diseases and vary with farming practices in the Mediterranean Thau lagoon, France. Aquacult Environ Interact 2:215-237

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

Pernet F, Gachelin S, Stanisière JY, Petton B, Fleury E, Mazurié J, Byron C (2019) Farmer monitoring reveals the effect of tidal height on mortality risk of oysters during a herpesvirus outbreak. ICES J Mar Sci 76:1816-1824

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

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

R Core Team (2020) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna

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

Ren JS, Ross AH (2001) A dynamic energy budget model of the Pacific oyster *Crassostrea gigas*. Ecol Modell 142: 105-120

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

Ren JS, Ross AH, Hadfield MG, Hayden BJ (2010) An ecosystem model for estimating potential shellfish culture production in sheltered coastal waters. Ecol Modell 221:527-539

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

Rosland R, Strand Ø, Alunno-Bruscia M, Bacher C, Strohmeier T (2009) Applying Dynamic Energy Budget (DEB) theory to simulate growth and bio-energetics of blue mussels under low seston conditions. J Sea Res 62: 49-61

Saulquin B, Gohin F, Garrello R (2011) Regional objective analysis for merging high-resolution MERIS, MODIS/Aqua, and SeaWiFS chlorophyll-a data from 1998 to 2008 on the European Atlantic shelf. IEEE Trans Geosci Rem Sens 49:143-154

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

Smith EL (1936) Photosynthesis in relation to light and carbon dioxide. Proc Natl Acad Sci USA 22:504-511
<figure><img src="image_70.png" /></figure>

Steeves LE, Filgueira R, Guyondet T, Chassé J, Comeau L (2018) Past, present, and future: performance of two bivalve species under changing environmental conditions. *Front Mar Sci* 5:184

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

Thomas Y, Bacher C (2018) Assessing the sensitivity of bivalve populations to global warming using an individual-based modelling approach. *Glob Change Biol* 24: 4581–4597

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

Thomas Y, Mazurié J, Alunno-Bruscia M, Bacher C and others (2011) Modelling spatio-temporal variability of *Mytilus edulis* (L.) growth by forcing a dynamic energy budget model with satellite-derived environmental data. *J Sea Res* 66:308–317

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

Thomas Y, Pouvreau S, Alunno-Bruscia M, Barillé L, Gohin F, Bryère P, Gernez P (2016) Global change and climate-driven invasion of the Pacific oyster (*Crassostrea gigas*) along European coasts: a bioenergetics modelling approach. *J Biogeogr* 43:568–579

Thomas Y, Cassou C, Gernez P, Pouvreau S (2018) Oysters as sentinels of climate variability and climate change in coastal ecosystems. *Environ Res Lett* 13:104009

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

Thrush MA, Pearce FM, Gubbins MJ, Oldtmann BC, Peeler EJ (2017) A simple model to rank shellfish farming areas based on the risk of disease introduction and spread. *Transbound Emerg Dis* 64:1200–1209

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

Tissot C, Brosset D, Barillé L, Le Grel L, Tillier I, Rouan M, Le Tixerant M (2012) Modeling oyster farming activities in coastal areas: a generic framework and preliminary application to a case study. *Coast Manage* 40:484–500

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

Vanhoutte-Brunier A, Fernand L, Ménesguen A, Lyons S, Gohin F, Cugier P (2008) Modelling the *Karenia mikimotoi* bloom that occurred in the western English Channel during summer 2003. *Ecol Modell* 210:351–376

Weitzman J, Filgueira R (2020) The evolution and application of carrying capacity in aquaculture: towards a research agenda. *Rev Aquacult* 12:1297–1322

Editorial responsibility: Gianluca Sarà,
Palermo, Italy

Reviewed by: D. Brigolin, R. Filgueira

Submitted: June 17, 2021

Accepted: January 27, 2022

Proofs received from author(s): April 1, 2022