# Should we account for mesozooplankton reproduction and ontogenetic growth in biogeochemical modeling?

Corentin Clerc, Olivier Aumont, Laurent Bopp

► To cite this version:

Corentin Clerc, Olivier Aumont, Laurent Bopp. Should we account for mesozooplankton reproduction and ontogenetic growth in biogeochemical modeling?. Theoretical Ecology, 2021, 10.1007/s12080-021-00519-5. [https://doi.org/10.1007/s12080-021-00519-5](https://doi.org/10.1007/s12080-021-00519-5)

HAL Id: hal-03318248

https://hal.science/hal-03318248

Submitted on 10 Aug 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
# Should we account for mesozooplankton reproduction and ontogenetic growth in biogeochemical modeling?

Corentin Clerc<fnref n="1" /> ID · Olivier Aumont<fnref n="2" /> ID · Laurent Bopp<fnref n="1" /> ID

Received: 28 May 2020 / Accepted: 11 June 2021

© The Author(s) 2021

## Abstract

Mesozooplankton play a key role in marine ecosystems as they modulate the transfer of energy from phytoplankton to large marine organisms. In addition, they directly influence the oceanic cycles of carbon and nutrients through vertical migrations, fecal pellet production, respiration, and excretion. Mesozooplankton are mainly made up of metazoans, which undergo important size changes during their life cycle, resulting in significant variations in metabolic rates. However, most marine biogeochemical models represent mesozooplankton as protists-like organisms. Here, we study the potential caveats of this simplistic representation by using a chemostat-like zero-dimensional model with four different Nutrient-Phytoplankton-Zooplankton configurations in which the description of mesozooplankton ranges from protist-type organisms to using a size-based formulation including explicit reproduction and ontogenetic growth. We show that the size-based formulation strongly impacts mesozooplankton. First, it generates a delay of a few months in the response to an increase in food availability. Second, the increase in mesozooplankton biomass displays much larger temporal variations, in the form of successive cohorts, because of the dependency of the ingestion rate to body size. However, the size-based formulation does not affect smaller plankton or nutrient concentrations. A proper assessment of these top-down effects would require implementing our size-resolved approach in a 3-dimensional biogeochemical model. Furthermore, the bottom-up effects on higher trophic levels resulting from the significant changes in the temporal dynamics of mesozooplankton could be estimated in an end-to-end model coupling low and high trophic levels.

**Keywords** Mesozooplankton · Biogeochemical models · Size-spectrum · Trait-based modeling · Reproduction

## Introduction

Zooplankton plays a pivotal role in the marine ecosystems because they are at the interface between primary producers and higher trophic levels. They are as such responsible for a large fraction of the energy transfer to fish (Beaugrand et al. 2010; Carlotti and Poggiale 2010; Mitra et al. 2014). They are also key players in marine biogeochemical cycles as they convert a large fraction of the organic matter they consume to dissolved inorganic carbon and nutrient pools. Furthermore, they modulate the amount of organic matter

that is exported to the interior of the ocean either through the sinking of fecal pellets and dead organisms (Henschke et al. 2016; Steinberg and Landry 2017) and through diel and seasonal vertical migration (Emerson et al. 1997; Packard and Gomez 2013; Aumont et al. 2018).

Assessment of trophic flows through zooplankton is difficult to perform due to insufficient data and methodological difficulties (Buitenhuis et al. 2006; Everett et al. 2017). Nevertheless, several estimates on the global scale have been proposed for the magnitude of feeding by different zooplankton size classes on primary autotrophic producers, i.e., phytoplankton. Microzooplankton feeding is found to consume from 50% up to almost 80% of primary production (PP) and exerts thus a major top-down control on phytoplankton (Calbet and Landry 2004; Schmoker et al. 2013). Mesozooplankton respiration has been estimated as 10% to about 30% of global PP, both phytoplankton and microzooplankton contributing importantly to the diet of mesozooplankton (Calbet 2001; Hernandez-Len and Ikeda 2005).

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

<footnote>Corentin Clerc
corentin.clerc@lmd.ens.fr</footnote>

<footnote n="1">¹ LMD / IPSL, Ecole normale supérieure / Université PSL, CNRS, Ecole Polytechnique, Sorbonne Université, Paris, France</footnote>

<footnote n="2">² LOCEAN / IPSL, IRD, CNRS, Sorbonne Université, MNHN, Paris, France</footnote>
A first separation that should be made within zooplankton is the distinction between protozoan unicellular organisms which can be referred to as microzooplankton and metazoan pluricellular organisms which are in general larger animals included in mesozooplankton and macrozooplankton (Carlotti and Poggiale 2010; Mitra et al. 2014). Despite the fact that protozoan zooplankton can have quite complex lifecycles, they share many similarities with phytoplankton except for their heterotrophy. Their size changes are limited to about doubling or halving of their cell volume, which involves small changes in their metabolic rates during their life history. On the other hand, metazoan zooplankton undergo large changes in their size, up to several orders of magnitude, over their life cycle. As a consequence, their metabolic rates change strongly over their life cycle history. Furthermore, intraguild predation (cannibalism) is common in mesozooplankton, for instance among copepods and chaetognaths (Brodeur and Terazaki 1999; Basedow and Tande 2006; Ohman et al. 2008).

Mathematical models of marine biogeochemistry and
plankton dynamics have gained a considerable interest over
the last decades and are now a central plank in ecological,
biogeochemical, and climate change studies. Most of the
biogeochemical models used today find their roots in the
NPZ formalism popularized by Fasham et al. (1990). Their
complexity has strongly increased by adding multiple lim-
iting nutrients and multiple phytoplankton and zooplankton
functional groups or size classes, e.g., Moore et al. (2002);
Le Queré et al. (2005); Follows et al. (2007); Ward et al.
(2012); Aumont et al. (2015). However, most of the effort has
focused on the description of the nutrient-phytoplankton link.
In comparison, the representation of zooplankton has received
much less attention and is very often restricted to employing a
few size-classes, generally only 2 corresponding to micro- and
mesozooplankton, both considered as protist-like organisms
(Sailley et al. 2013). Indeed, for the vast majority of past and
current biogeochemical models, zooplankton act as a closure
term for nutrient and carbon fluxes.

Physiologically Structured Population Models (PSPMs) (Metz and Diekmann 1986; De Roos 1997) is a class of models that are designed to investigate the dynamics of population controlled by size-dependent processes. PSPMs are based on a description of relevant processes at the individual level such as feeding, development, reproduction and mortality which depend on the environment (for instance, the amount of food) and on the state of the organism itself. As size controls most, if not all, of these processes, including feeding interactions, PSMPs are often constructed as size-resolved models, (e.g., Persson et al. 1998; van Kooten et al. 2007). Indeed, this subclass of PSPMs has been used to improve our understanding of the effects of ontogenetic growth at the population level (De Roos 2018) as well as at the community level (De Roos 2020). In systems where size-structured predators feed on unstructured preys, cyclic cohort cycles may appear with periods

depending on the competitive advantage between juveniles
and adults (De Roos et al. 2003; De Roos and Persson 2003;
De Roos et al. 2008). This type of model may also exhibit
alternative stable states, at least when adult and juvenile stages
feed on different resources (Guill 2009).

The complexity of size-structured PSPMs makes their
application at the community level challenging, because they
require a large set of species-specific parameters. As a more
tractable alternative, aggregated stage-based approaches have
been derived from PSPMs (De Roos et al. 2008; Guill 2009;
Soudijn and De Roos 2017). However, even if they proved to be
able to capture some of the properties of the PSPMs, this is not
always the case and their predictions depend on the number of
stages they represent (Guill 2009). Purely size-based models
have been proposed to represent communities which reduce
the complexity of PSPMs by assigning a single trait to dif-
ferentiate between species. This trait can be size at maturation
(Hartvig et al. 2011) or maximum size (Maury and Poggiale
2013). The structuring role of size at the community level is
supported by marine observations that often show continuous
and quasi-linear size spectra in logarithmic space (Sheldon
et al. 1972; Quinones et al. 2003; Zhou et al. 2009; Thompson
et al. 2013). The dynamics of the abundance distribution is
generally described by the well-known McKendrick von
Foerster equation that combines an advection term along the
size dimension corresponding to growth and some sink terms
due to predation and other mortality processes (Guiet et al.
2016). Using these different approaches, numerous models at
the population or community levels have been developed to
study zooplankton and their response to both biotic and abiotic
changes in their environment (e.g., Carlotti et al. 2000; Flynn
and Irigoien 2009; Zhou et al. 2010; Pinceel et al. 2016). Some
of these models have been integrated in biogeochemical mod-
els (e.g., Fennel 2001; Baird and Suthers 2007; Zhou et al.
2010; Serra-Pompei et al. 2020). However, to our knowledge,
the consequences of explicitly describing ontogenetic growth
and reproduction of mesozooplankton have not been exten-
sively explored, especially concerning biogeochemistry and
plankton dynamics. Furthermore, no spatially resolved biogeo-
chemical modeling study includes such a detailed description
of mesozooplankton community.

In this paper, we propose to investigate the consequences of
an explicit description of individual growth and reproduction of
a mesozooplankton community on the temporal dynamics of a
simple NPZ-type plankton model. To achieve this goal, we have
developed a suite of four different model configurations in which
the description of mesozooplankton goes from the dominant rep-
resentation in biogeochemical models, as protist-like organisms, to
a consistent size-based formulation with reproduction and ontoge-
netic growth. These different model configurations are used in a
zero-dimensional chemostat setting in which the incoming nutri-
ent levels are varied with time in order to assess the impact on the
temporal dynamics of nutrient and plankton biomass.

<figure><img src="image_3.png" /></figure>
# Materials and method

## Model description

**Model structure** This study is based on the simple biogeochemical model developed by Mayersohn (2018). The state variables of this model represent a resource or a group of organisms whose temporal evolution is governed by an ordinary differential equation (ODE). These state variables are expressed in terms of their phosphorus content (µmolP L⁻¹). The original version is a classical NPZ model with three major types of compartments: nutrient, phytoplankton, and zooplankton. The version developed for this study includes a subdivision of the zooplankton compartment to explicitly represent mesozooplankton size classes, and thus sexual reproduction and growth. Figure 1 summarizes the different (sub)compartments of the model and their interactions.

The model contains 4 main compartments (dotted round squares in Fig. 1), expressed in terms of biomass or concentration, and each comprising one or more subcompartments (simple round squares in Fig. 1). Phosphate is the only resource considered here. Phytoplankton *P* is divided into three size classes *P*ᵢ where *i* ∈ [0, 2]. Zooplankton is divided into two main compartments according to size: microzooplankton *Z* and mesozooplankton

*M*. Microzooplankton is further divided into strictly heterotrophic protists *U* (i.e., strict microzooplankton) and Ns/2 size classes representing juvenile multi-cellular zooplankton *J*ᵢ with *i* ∈ [0, Ns/2 - 1]. The mesozooplankton compartment is divided into Ns/2 mature life stages *A*ᵢ with *i* ∈ [0, Ns/2 - 1]. The maximum size mesozooplankton compartment is also referred to as *A*max. The symbols are summarized in Table 1.

Our system represents a chemostat experiment in a zero-dimensional (0-D) setting. This system can be considered similar to an idealized two-layer ocean model with instantaneous remineralization of organic matter under the mixed layer (Evans et al. 1985; Thingstad et al. 1996; Ward et al. 2012). Inorganic phosphorus (phosphate) is supplied to the system at a specified dilution rate τ(t) and concentration C(t), which can be either constant or vary with time, to mimic vertical mixing events with nutrient-rich deep waters. An equal outflow dilutes inorganic phosphorus and the organic living compartments in the system. Available Photosynthetically Active Radiation (PAR) *I(t)* is specified as an external environmental variable and is supposed homogeneous within the system. The external environmental conditions are represented by bold arrows in Fig. 1.

**General system dynamics** The following set of equations governs the temporal evolution of the inorganic resource and phytoplankton groups:

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Model structure. The dotted round squares are the major functional compartments. Plain round squares are the sub-compartments (i.e., Resources, Phytoplankton and Zooplankton). The thin arrows represent the flows between compartments or sub-compartments. Bold arrows represent the external forcings. The large round white dotted square is the modeled system, which is similar to a chemostat. The small round blue dotted square is the external medium with which our system is mixed. “Repro.” claims for “Reproduction”</figcaption></figure>
$$ \frac{dR}{dt} = \tau(C(t) - R) - \delta \underbrace{\sum_{i=0}^{2} r_{P_i} P_i}_{\text{mixing}} + \underbrace{mM^2}_{\text{mortality}} + \underbrace{(1 - \beta)(r_M M + r_Z Z)}_{\text{sloppy feeding}} + \underbrace{(1 - w)\beta r_{A_{\text{max}}} A_{\text{max}}}_{\text{growth closure}} $$

$$ \frac{dP_i}{dt} = \tau(Q_{\text{min}} - P_i) + \underbrace{\left[ \delta r_{P_i} \underbrace{- \rho_{P_i, Z} \frac{1}{T_Z} r_Z Z + \rho_{P_i, M} \frac{1}{T_M} r_M M}_{\text{predation}} \right]}_{\text{mixing}} \cdot P_i $$

Growth rate of phytoplankton $r_{P_i}$ depends on light and on nutrient availability (Table 1). The only loss term for phytoplankton is predation by the different zooplankton groups. The inorganic resource $R$ is taken up by phytoplankton and is resupplied by mortality and excretion from zooplankton.

The temporal evolution of the different zooplankton groups is computed according to the following set of equations:

The formulation for the different terms used in this set of equations is detailed in Table 1 (Parts IV and V). All zooplankton groups feed on all three phytoplankton size classes. In addition, mesozooplankton feed on strict microzooplankton and on juveniles. Grazing is a function of the amount of food available ($T_Z$ for strict microzooplankton and juveniles, $T_M$ for mesozooplankton). It is parameterized according to a Holling type III function that improves the stability of the model (Gentleman and Neuheimer 2008). Feeding on

$$ \frac{dU}{dt} = \tau(Q_{\text{min}} - U) + \underbrace{\left[ \beta r_U - \frac{1}{T_M} r_M M \right]}_{\text{growth}} \cdot U $$

$$ \frac{dJ_0}{dt} = \underbrace{-\tau}_{\text{mixing}} + \underbrace{(1 - \nu)\beta r_{J_0}}_{\text{growth and transition}} - \underbrace{\frac{1}{T_M} r_M M}_{\text{predation}} \cdot J_0 + \underbrace{w\beta r_{A_0} A_0}_{\text{reproduction}} $$

$$ \frac{dJ_s}{dt} = \underbrace{-\tau}_{\text{mixing}} + \underbrace{(1 - \nu)\beta r_{J_s}}_{\text{growth and transition}} - \underbrace{\frac{1}{T_M} r_M M}_{\text{predation}} \cdot J_s + \underbrace{\nu\beta r_{J_{s-1}} J_{s-1}}_{\text{transition}} + \underbrace{w\beta r_{A_s} A_s}_{\text{reproduction}} $$

$$ \frac{dA_0}{dt} = \underbrace{-\tau}_{\text{mixing}} + \underbrace{(1 - w)(1 - \nu)\beta r_{A_0}}_{\text{growth, reproduction and transition}} - \underbrace{mM}_{\text{mortality}} \cdot A_0 + \underbrace{\nu\beta r_{J_{N_s/2-1}} J_{N_s/2-1}}_{\text{transition}} $$

$$ \frac{dA_s}{dt} = \underbrace{-\tau}_{\text{mixing}} + \underbrace{(1 - w)(1 - \nu)\beta r_{A_s}}_{\text{growth, reproduction and transition}} - \underbrace{mM}_{\text{mortality}} \cdot A_s + \underbrace{(1 - w)\nu\beta r_{A_{s-1}} A_{s-1}}_{\text{transition}} $$

$$ \frac{dA_{\text{max}}}{dt} = \tau(Q_{\text{min}} - A_{\text{max}}) - \underbrace{mM A_{\text{max}}}_{\text{mortality}} + \underbrace{(1 - w)\nu\beta r_{A_{N_s/2-2}} A_{N_s/2-2}}_{\text{transition}} $$
<table><caption>**Table 1** Variables and parameters used in the set of equations governing the temporal evolution of the state variables</caption><thead><tr><th>Symbol</th><th>Expression (if applicable)</th><th>Description</th></tr></thead><tbody><tr><td colspan="3"><strong>I. STATE VARIABLES</strong></td></tr><tr><td>R</td><td></td><td>Inorganic resource (here PO<sub>4</sub><sup>2-</sup>)</td></tr><tr><td>P<sub>i</sub></td><td></td><td>Phytoplankton size-class i</td></tr><tr><td>U</td><td></td><td>Unicellular microzooplankton</td></tr><tr><td>J<sub>s</sub></td><td></td><td>Juvenile size-class s of mesozooplankton</td></tr><tr><td>A<sub>s</sub></td><td></td><td>Mature size-class s of mesozooplankton</td></tr><tr><td>A<sub>max</sub> = A<sub>N<sub>s</sub>/2-1</sub></td><td></td><td>Largest mesozooplankton size-class</td></tr><tr><td>N<sub>s</sub></td><td></td><td>Number of mesozooplankton size-classes</td></tr><tr><td>P</td><td>= &sum;<sub>i=0</sub><sup>2</sup> P<sub>i</sub></td><td>Total biomass of phytoplankton</td></tr><tr><td>Z</td><td>= U + &sum;<sub>s=0</sub><sup>N<sub>s</sub>/2-1</sup> J<sub>s</sub></td><td>Total biomass of microzooplankton</td></tr><tr><td>M</td><td>= &sum;<sub>s=0</sub><sup>N<sub>s</sub>/2-1</sup> A<sub>s</sub></td><td>Total biomass of mesozooplankton</td></tr><tr><td colspan="3"><strong>II. EXTERNAL FORCINGS</strong></td></tr><tr><td>&tau;</td><td></td><td>Dilution rate</td></tr><tr><td>T</td><td></td><td>Period of the chemostat step function</td></tr><tr><td>C<sub>0</sub></td><td></td><td>Maximum resource concentration in the inflow</td></tr><tr><td>&theta;</td><td></td><td>Nutrient to carbon ratio in organic matter</td></tr><tr><td>C(t)</td><td>= C<sub>0</sub> &int;<sub>{t &isin; [kT, kT + &tau;/2, 1, k &isin; N]}</sub> &int;<sub>{t &gt; 1yr}</sub> 1</td><td>Resource concentration in the inflow</td></tr><tr><td>I</td><td></td><td>Light intensity</td></tr><tr><td>Q<sub>min</sub></td><td></td><td>Organisms threshold</td></tr><tr><td colspan="3"><strong>III. PHOTOSYNTHESIS</strong></td></tr><tr><td>K<sup>L</sup></td><td></td><td>Light saturation constant of phytoplankton</td></tr><tr><td>&psi;</td><td>= 1 - e<sup>-I/kL</sup></td><td>Light limitation factor of phytoplankton</td></tr><tr><td>K<sup>N</sup><sub>P<sub>i</sub></sub></td><td></td><td>Nutrient half-saturation constant of phytoplankton P<sub>i</sub></td></tr><tr><td>S<sub>P<sub>i</sub></sub></td><td>= R / (K<sup>N</sup><sub>P<sub>i</sub></sub> + R)</td><td>Nutrient limitation factor of phytoplankton P<sub>i</sub></td></tr><tr><td>g<sub>P</sub></td><td></td><td>Maximum growth rate of phytoplankton</td></tr><tr><td>r<sub>P<sub>i</sub></sub></td><td>= g<sub>P</sub>&psi;S<sub>P<sub>i</sub></sub></td><td>Growth rate of phytoplankton P<sub>i</sub></td></tr><tr><td>&delta;</td><td></td><td>Growth efficiency</td></tr><tr><td colspan="3"><strong>IV. PREDATION</strong></td></tr><tr><td>&rho;<sub>P<sub>i</sub>,Z</sub></td><td></td><td>Microzooplankton preference for phytoplankton P<sub>i</sub></td></tr><tr><td>&rho;<sub>P<sub>i</sub>,M</sub></td><td></td><td>Mesozooplankton preference for phytoplankton P<sub>i</sub></td></tr><tr><td>T<sub>Z</sub></td><td>= &sum;<sub>i=0</sub><sup>2</sup> &rho;<sub>P<sub>i</sub>,Z</sub>P<sub>i</sub></td><td>Amount of food available for microzooplankton</td></tr><tr><td>T<sub>M</sub></td><td>= &sum;<sub>i=0</sub><sup>2</sup> &rho;<sub>P<sub>i</sub>,M</sub>P<sub>i</sub> + Z</td><td>Amount of food available for mesozooplankton</td></tr><tr><td>K<sup>G</sup></td><td></td><td>Food half-saturation constant of zooplankton</td></tr><tr><td>S<sub>Z</sub></td><td>= T<sub>Z</sub> / (K<sup>G2</sup>+T<sub>Z</sub><sup>2</sup>)</td><td>Food functional response of microzooplankton</td></tr><tr><td>S<sub>M</sub></td><td>= T<sub>M</sub> / (K<sup>G2</sup>+T<sub>M</sub><sup>2</sup>)</td><td>Food functional response of mesozooplankton</td></tr><tr><td>g<sub>x</sub></td><td></td><td>Maximum ingestion rate of zooplankton x &isin; Z &cup; M</td></tr><tr><td>g<sup>e</sup><sub>Z</sub></td><td>= 1/Z &sum;<sub>z &isin; Z</sub> zg<sub>z</sub></td><td>Mean maximum ingestion rate of microzooplankton</td></tr><tr><td>g<sup>e</sup><sub>M</sub></td><td>= 1/M &sum;<sub>m &isin; M</sub> mg<sub>m</sub></td><td>Mean maximum ingestion rate of mesozooplankton</td></tr><tr><td>r<sub>Z</sub></td><td>= g<sub>Z</sub>S<sub>Z</sub></td><td>Ingestion rate of microzooplankton compartment z &isin; Z</td></tr></tbody></table>
<table><caption>**Table 1** (continued)</caption><thead><tr><th>Symbol</th><th>Expression (if applicable)</th><th>Description</th></tr></thead><tbody><tr><td>r<sub>m</sub></td><td>= g<sub>m</sub>S<sub>M</sub></td><td>Ingestion rate of mesozooplankton compartment m ∈ M</td></tr><tr><td>r<sub>Z</sub></td><td>= g<sub>Z</sub>S<sub>Z</sub></td><td>Mean ingestion rate of microzooplankton</td></tr><tr><td>r<sub>M</sub></td><td>= g<sub>M</sub><sup>e</sup>S<sub>M</sub></td><td>Mean ingestion rate of mesozooplankton</td></tr><tr><td>β</td><td></td><td>Growth efficiency of zooplankton</td></tr><tr><td>m</td><td></td><td>Quadratic mortality rate of mature mesozooplankton size-classes<br>V. REPRODUCTION AND GROWTH</td></tr><tr><td>w</td><td></td><td>Fraction of the assimilated food allocated to reproduction</td></tr><tr><td>v</td><td></td><td>Transition rate across metazoan size-classes</td></tr></tbody></table>

the different available preys is described by a proportional formulation which implies no-switching (Gentleman et al. 2003).

For each mature mesozooplankton $A_s$, part of the assimilated food $w$ is allocated to reproduction and is transferred to the juvenile sub-compartment $J_s$. The remainder of the assimilated food is used for growth, resulting in a transfer between adjacent size classes at a rate $v$. The value of this parameter depends on the number of size classes and the assumed size distribution within each size class (see Supplementary Material Section 1). For the largest size class of mature mesozooplankton $A_{max}$, no size growth is possible. Thus, the amount of assimilated food that would lead to an increase in size is channeled to the inorganic resource (we refer to this term as growth closure). Furthermore, we assume that the deepwater mass contains only mature life stages for mesozooplankton.

The overall dynamics of mesozooplankton $M$ (Eq. 2) and microzooplankton $Z$ (Eq. 1) compartments, on which we focus in this study, is obtained by summing the terms of the different sub-compartments:

$$ \frac{dZ}{dt} = \tau(Q_{min} - Z) + \underbrace{\beta r_Z Z}_{\text{net growth}} - \underbrace{\frac{Z}{T_M} r_M M}_{\text{predation}} + \underbrace{w \beta r_M M}_{\text{reproduction}} - \underbrace{v \beta r_{J_{N_s/2-1}} J_{N_s/2-1}}_{\text{maturation}} \quad (1) $$

$$ \frac{dM}{dt} = \tau(Q_{min} - M) + \underbrace{(1 - w)\beta r_M M}_{\text{net growth}} - \underbrace{mM^2}_{\text{mortality}} + \underbrace{v\beta r_{J_{N_s/2-1}} J_{N_s/2-1}}_{\text{maturation}} - \underbrace{(1 - w)\beta r_{A_{max}} A_{max}}_{\text{growth closure}} \quad (2) $$

**Size-based parameterization** Phytoplankton is divided into three classes of increasing size. The nutrient

half-saturation constant increases with size (e.g., Eppley et al. 1969; Edwards et al. 2012). Maximum growth rate is kept constant with size, despite observational studies that suggest a unimodal size scaling of maximum phytoplankton growth rate (e.g., Wirtz 2011; Maran et al. 2013). Zooplankton preference for phytoplankton is also set to vary with size: microzooplankton groups are assumed to preferentially feed on small phytoplankton, while mesozooplankton compartments are more likely to feed on large phytoplankton. Food preference is constant in each major zooplankton compartment (microzooplankton and mesozooplankton).

We consider a subdivision of the two metazoan zooplankton compartments into size classes of equal width in logarithmic space. The center of each size class is defined as follows: $l_s = l_{min} \left( \frac{l_{max}}{l_{min}} \right)^{\frac{2s+1}{2N_s}}$ where $s \in [0, N_s - 1]$. The width of each size class is $\Delta \ln(l_s) = \frac{1}{N_s} \ln \left( \frac{l_{max}}{l_{min}} \right)$ in logarithmic space and is therefore constant. We define the length factor of the size class $X_s$ as $\frac{2s+1}{2N_s}$ which goes from 0 to 1.

The maximum ingestion rate of the different zooplankton classes is set according to an allometric relationship proposed by Hansen et al. (1997). The half-saturation constant used in the grazing parameterization is supposed constant as observations suggest no significant variations with size (Hansen et al. 1997). The transition rate $v$ between the different size-classes was computed by assuming that the slope of the biomass size-spectrum within each size-class is constant in a log-log space. It is set to -3 following the seminal study of Sheldon et al. (1972), which corresponds to an approximate constant biomass in logarithmically equal size intervals. The expressions for the transition rate and for the maximum ingestion rate are shown in Table 2 and their computation is detailed in the supplementary material (Sections 1 and 2). The size-dependent formulations used in our standard model configuration are listed in Table 2.
<table><caption>**Table 2** Parameters and equations used in the size-based parameterizations</caption><thead><tr><th>Term</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td></td><td></td><td>PHYTOPLANKTON</td></tr><tr><td>L(P<sub>i</sub>)</td><td>= <sup>i+1</sup>&frasl;<sub>4</sub></td><td>Length factor for phytoplankton P<sub>i</sub></td></tr><tr><td>K<sup>N</sup><sub>P<sub>i</sub></sub></td><td>= (K<sup>N</sup><sub>min</sub> + (K<sup>N</sup><sub>max</sub> - K<sup>N</sup><sub>min</sub>)L(P<sub>i</sub>)<sup>2</sup>)</td><td>Nutrient half-saturation constant for phytoplankton P<sub>i</sub></td></tr><tr><td></td><td></td><td>ZOOPLANKTON</td></tr><tr><td>l<sub>min</sub></td><td></td><td>Minimal metazoan zooplankton body length</td></tr><tr><td>l<sub>max</sub></td><td></td><td>Maximal metazoan zooplankton body length</td></tr><tr><td>v</td><td>= <sup>N<sub>s</sub></sup>&frasl;<sub>3 ln <sup>l<sub>max</sub></sup>&frasl;<sub>l<sub>min</sub></sub></sub></td><td>Transition rate between the mesozooplankton size-classes</td></tr><tr><td>ρ<sub>P<sub>i</sub>,Z</sub></td><td>= (1 - L(P<sub>i</sub>)<sup>2</sup>)</td><td>Microzooplankton preference for phytoplankton P<sub>i</sub></td></tr><tr><td>ρ<sub>P<sub>i</sub>,M</sub></td><td>= L(P<sub>i</sub>)<sup>2</sup></td><td>Mesozooplankton preference for phytoplankton P<sub>i</sub></td></tr><tr><td>g<sub>M</sub></td><td></td><td>Geometric mean of the maximum mesozooplankton ingestion rate</td></tr><tr><td>g<sub>Z</sub></td><td></td><td>Geometric mean of the maximum microzooplankton ingestion rate</td></tr><tr><td>L(J<sub>s</sub>)</td><td>= <sup>2s+1</sup>&frasl;<sub>2N<sub>s</sub></sub></td><td>Length factor of juvenile size-classes J<sub>s</sub></td></tr><tr><td>L(A<sub>s</sub>)</td><td>= <sup>N<sub>s</sub>+2s+1</sup>&frasl;<sub>2N<sub>s</sub></sub></td><td>Length factor of mature size-classes A<sub>s</sub></td></tr><tr><td>L(U)</td><td>= <sup>1</sup>&frasl;<sub>4</sub></td><td>Length factor for generic microzooplankton U</td></tr><tr><td>ln g<sub>s</sub></td><td>= ln g<sub>Z</sub> + α(L(U) - L(X<sub>s</sub>)) ln (<sup>l<sub>max</sub></sup>&frasl;<sub>l<sub>min</sub></sub>)</td><td>ALLOMETRIC SCALING OF RATES</td></tr><tr><td></td><td></td><td>Maximum ingestion rate of the zooplankton size-class X<sub>s</sub></td></tr><tr><td></td><td></td><td>CONSTANT RATES</td></tr><tr><td>g<sub>s</sub></td><td>= 1<sub>{X<sub>s</sub>∈Z}</sub>g<sub>Z</sub> + 1<sub>{X<sub>s</sub>∈M}</sub>g<sub>M</sub></td><td>Maximum grazing rate of the zooplankton size-class X<sub>s</sub></td></tr></tbody></table>

## Alternative models

The model described above is hereafter referred to as the “standard model”. It relies on an explicit description of the size distribution of metazoan zooplankton, including ontogenetic growth and reproduction. Furthermore, maximum ingestion rates vary with size following the allometric relationship proposed by Hansen et al. (1997). To investigate the effects of this representation of mesozooplankton ontogenetic growth and reproduction on the model behavior, three alternative model formulations were developed for comparison. The specificities of these models are presented in Fig. 2. Detailed equations are available in the supplementary materials (Section 3).

The first alternative model is a “constant-rate” model, in which maximum ingestion rates are constant across all size classes of each major zooplankton compartment. The geometric mean of the maximum ingestion rates of microzooplankton and mesozooplankton, respectively, is assigned to the juvenile and mature size classes, i.e., $g_Z$ and $g_M$ (Table 2). By comparison with the standard model, this model configuration is designed to evaluate the effects of allometric scaling of maximum ingestion rates on the temporal dynamics of the system.

The second alternative model is a “two-life-stage” model (Fig. 2b) in which there is no size discretization of juveniles and mature organisms. Consequently, the representation of metazoan zooplankton is limited to two life stages: juveniles

and mature organisms. In this model, adult mesozooplankton allocates part of their assimilated food to reproduction, which is represented as a flux of organic matter from mature mesozooplankton to the juvenile life stage. As a result of ontogenetic growth, juveniles are transferred to adult mesozooplankton at a rate $v$ computed according to the formula presented in Table 1. The value of $v$ is here less than 1 in contrast to the first two model configurations that have a better resolved size distribution. Thus, food assimilation leads to an increase in the biomass of each metazoan compartments. As mature mesozooplankton is represented by a single compartment, the specific treatment of the largest size-class as applied in the standard and constant-rate models is not used here. Comparing the standard model with the two-life-stage model helps to investigate the effect of an increased size resolution on the temporal dynamics of the system.

The last alternative model is a “no-life-stage” model (Fig. 2a), in which metazoans are represented by a single compartment corresponding to mesozooplankton. Thus juvenile and mature organisms are supposed to have the same metabolic rates and the same predation behavior. There is no internal predation within the mesozooplankton compartment. This model is very similar to the classical approach adopted by most of the PFT-type biogeochemical models. In this case, the representation of both microzooplankton and mesozooplankton is similar and corresponds to a formalism used for protists whose typical mode of reproduction is based on cell division.
a. No life-stage model

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

b. Two life-stages model

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

C. Constant vs allometric rate models

<figure><img src="image_7.png" /><figcaption>**Fig. 2** Standard and alternative model configurations used in this study. **a.** “no-life-stage” model configuration. **b.** “two-life-stage” model configuration. Left red arrow represents maturation, right red arrow represents reproduction. In a. and b., light grey items denote the unchanged parts of the model presented in Fig. 1. **c.** Parameterization</figcaption></figure>

tions of the maximum ingestion rate of mesozooplankton. Plain lines with dots represent the standard model parameterization (allometric), and dashed lines with squares represent a configuration in which the maximum ingestion rate is constant within each main compartment. Vertical grey lines are delimiting the 20 size-classes

**Numerical experiments**

Each model configuration was run for 4 years using a simple Euler forward integration scheme with a time step of 1 hour. Light intensity, organism concentrations in the incoming waters and the dilution rate are kept constant throughout the simulations. Light intensity is set to a non-limiting value of 60 W m⁻². Dilution rate is maintained at 0.1 d⁻¹. Organisms concentrations in the inflow waters are set to a very small value to avoid any extinction of communities. The temporal dynamics of the different model configurations is investigated by periodically varying the incoming resource concentration according to a step function with a period *T* of 2 years. Incoming nutrients are varied between a null value and a non-limiting concentration of 1 µmolP L⁻¹.

The model parameters and their values are listed in Table 3. In the standard and constant-rate models, metazoans are split into 20 size-classes, equally distributed between juveniles and mature organisms. The reproduction parameter is set to 0.3 (Kooijman 2013). The size parameterization defines metazoans in the model as a community of individuals ranging from 10 to 4000 µm with an mean egg to adult ratio of 1/20. All the other parameter values have been assigned to typical values found in the literature. Reference for each parameter is provided in Table 3. Some parameters of the standard model (denoted by a ★ in Table 3) have been varied by ±50% around their standard value to produce a series of sensitivity simulations discussed in the last section of the results.

<figure><img src="image_8.png" /></figure>
<table><caption>**Table 3** Parameter values of the default configuration. The stars indicate the parameters on which sensitivity analyses have been performed</caption><thead><tr><th>Parameter</th><th>Default</th><th>Unit</th><th>Description</th><th>Range</th><th>Source</th></tr></thead><tbody><tr><td>N<sub>S</sub></td><td>20</td><td>-</td><td>Number of mesozooplankton size-classes</td><td></td><td></td></tr><tr><td>τ</td><td>0.1</td><td>d<sup>-1</sup></td><td>Dilution rate</td><td></td><td></td></tr><tr><td>T</td><td>720</td><td>d</td><td>Period of the chemostat step function</td><td></td><td></td></tr><tr><td>C<sub>0</sub></td><td>1</td><td>μmolP L<sup>-1</sup></td><td>Maximum resource concentration in the inflow</td><td></td><td></td></tr><tr><td>I</td><td>60</td><td>W m<sup>-2</sup></td><td>Light intensity</td><td></td><td></td></tr><tr><td>Q<sub>min</sub></td><td>10<sup>-6</sup></td><td>μmolP L<sup>-1</sup></td><td>Organisms threshold</td><td></td><td></td></tr><tr><td>K<sup>L</sup></td><td>20</td><td>W m<sup>-2</sup></td><td>Light saturation constant of phytoplankton</td><td>2.5-75</td><td>(Edwards et al. 2016)</td></tr><tr><td>K<sup>N</sup><sub>min</sub></td><td>0.01</td><td>μmolP L<sup>-1</sup></td><td>Minimum nutrient half-saturation constant for phytoplankton</td><td>0.002-0.02</td><td>(Eppley et al. 1969; Chisholm 1992)</td></tr><tr><td>K<sup>N</sup><sub>max</sub></td><td>0.08</td><td>μmolP L<sup>-1</sup></td><td>Maximum nutrient half-saturation constant for phytoplankton</td><td>0.02-0.2</td><td>(Eppley et al. 1969; Chisholm 1992)</td></tr><tr><td>g<sub>P</sub></td><td>2.4</td><td>d<sup>-1</sup></td><td>Maximum growth rate of phytoplankton</td><td>0.75-3.0</td><td>(Eppley 1972; Bissinger et al. 2008)</td></tr><tr><td>δ</td><td>0.95</td><td>-</td><td>Growth efficiency of phytoplankton</td><td></td><td>(Aumont et al. 2015)</td></tr><tr><td>K<sup>G</sup></td><td>0.04</td><td>μmolP L<sup>-1</sup></td><td>Food half-saturation constant of zooplankton</td><td>0.03-0.15</td><td>(Hansen et al. 1997; Stock and Dunne 2010)</td></tr><tr><td>β</td><td>0.4</td><td>-</td><td>Growth efficiency of zooplankton</td><td>0.1-0.5</td><td>(Buitenhuis et al. 2006; Straile 1997)</td></tr><tr><td>g<sub>M</sub></td><td>0.47</td><td>d<sup>-1</sup></td><td>Geometric mean of the maximum mesozooplankton ingestion rate</td><td>0.13-0.97</td><td>(Buitenhuis et al. 2006)</td></tr><tr><td>g<sub>Z</sub></td><td>2.0</td><td>d<sup>-1</sup></td><td>Geometric mean of the maximum microzooplankton ingestion rate</td><td>0.55-4.1</td><td>See table 2</td></tr><tr><td>m</td><td>3.0</td><td>μmolP L<sup>-1</sup>d<sup>-1</sup></td><td>Quadratic mortality rate of mature mesozooplankton</td><td></td><td>(Aumont et al. 2015)</td></tr><tr><td>w</td><td>0.3</td><td>-</td><td>Fraction of the assimilated food allocated to reproduction</td><td>0.2-0.8</td><td>(Kooijman 2013)</td></tr><tr><td>v</td><td>1.1</td><td>-</td><td>Transition rate across metazoan size-classes</td><td></td><td>See table 2</td></tr><tr><td>l<sub>min</sub></td><td>10</td><td>μm</td><td>Minimal metazoan zooplankton body length</td><td></td><td></td></tr><tr><td>l<sub>max</sub></td><td>4000</td><td>μm</td><td>Maximal metazoan zooplankton body length</td><td></td><td></td></tr><tr><td>α</td><td>0.48</td><td>-</td><td>Allometric parameter</td><td>0.42-0.54</td><td>(Hansen et al. 1997)</td></tr></tbody></table>

An analytical study of the size spectra at equilibrium allows a comparison of the theoretical continuous spectra with the simulated discrete spectra. This comparison presented in the supplementary materials (Section 4) confirms that 1) the discretization into 20 size classes is sufficient enough to derive a good representation of the size growth dynamics according to the McKendrick Von Foerster equation and, 2) the formulation of the transition rate *v* is robust to a deviation from the assumption that size spectrum has a slope of -3 in log-log space.

## Results

For all model simulations, we only analyze the last two years of the 4-yr simulations. Hereafter, we refer to the period following the switch from No-Nutrient Mode to

High-Nutrient Mode as *HNM* and the period following the switch from High-Nutrient Mode to No-Nutrient Mode as *ONM*.

In the No-Nutrient mode (*ONM*), the temporal dynamics across all model versions is similar (Fig. 3), with a sharp decrease for all compartments (resource and biomass) towards low values after ~50 days. This response is dictated by the mixing term with inflow waters that contain no nutrients and virtually no phytoplankton and zooplankton. The delay in the response varies depending on the compartment, i.e., more rapid for nutrient concentrations and less so for phytoplankton and zooplankton biomasses.

Hereafter, we focus our analysis on the response to the High-Nutrient mode (*HNM*) and present the results of the standard model first, and then a comparison to the alternative models.
<figure><img src="image_9.png" /><figcaption>**Fig. 3** Time evolution of the phosphorus content of all simulated compartments for all model versions. White area represents the high-nutrient mode, orange areas the no-nutrient mode. Dashed lines in the microzooplankton panel represent the juvenile metazoans</figcaption></figure>

**Standard model**

The standard model output is represented in orange on Fig. 3. In *HNM*, about 70% (i.e. 0.7 µmol P L⁻¹) of the total phosphorus content of the system is not converted into biomass at equilibrium, so that phytoplankton growth is not limited by nutrient. Contrary to *ONM*, the temporal dynamics varies a lot across compartments, with equilibrium values

reached much earlier for nutrients, phytoplankton, and microzooplankton biomasses than for mesozooplankton.

A short phytoplankton bloom is observed during the first 10 days. During this phase, up to 65% of the total phosphorus (*P*) content of the system is converted into phytoplankton biomass. Following the phytoplankton bloom, microzooplankton biomass also increases sharply reaching 55% of the total *P* content by the end of the phytoplankton bloom.

<figure><img src="image_10.png" /></figure>
<figure><img src="image_11.png" /><figcaption>**Fig. 4** Log-time evolution of the relative biomass distribution within size-classes for the standard model. Red dots indicate mature metazoan size-classes (in the mesozooplankton compartment), pink dots indicate juvenile metazoan size-classes (in the microzooplankton compartment)</figcaption></figure>

The high grazing pressure exerted by microzooplankton leads to a rapid collapse of phytoplankton at $t = 10$ days. Microzooplankton biomass then decreases, being limited by prey abundance, and reaches an equilibrium ~ 1 month after the switch to the high-nutrient mode. To summarize, both microzooplankton and phytoplankton experience a rapid and intense early bloom phase that is then followed by a rapid decline. Apart from the mesozooplankton compartment, a quasi-equilibrium is reached by the end of the first month in *HNM*.

For mesozooplankton, the temporal dynamics is markedly different from that of phyto- or microzooplankton. No outburst is observed and the temporal evolution is indeed much slower (Fig. 3). Consequently, equilibrium takes several months (> 1 year) to be reached.

To get some insights into the mesozooplankton response to *HNM*, we explore the temporal dynamics of the 20 different size-classes of metazoans, for both juveniles (in the microzooplankton compartment) and mature individuals (in the mesozooplankton compartment). Figure 4 shows the evolution of the relative abundance of biomass in the different size-classes over time (i.e., the size-spectrum time evolution). During the first 9 days of *HNM*, biomass is concentrated in the mesozooplankton size classes, as the inflow waters contain only mature organisms.

Reproduction of these mature organisms produces juveniles and induces the propagation of a cohort across the size spectrum as individuals keep progressing in their size-classes. This first cohort reaches the maximum size-class after about two months in *HNM*. When this cohort enters the mature size-classes domain (red points on Fig. 4), it triggers reproduction that leads to the production of new juveniles and thus, generates a second cohort. During that second cohort, the relative biomass of the mature size-classes reaches lower values than the first cohort because of the progressive build-up of juveniles (Fig. 3). The start of a third and a fourth cohort is observed before the system reaches equilibrium.

The temporal evolution of mesozooplankton in *HNM* can be split into 5 phases (Fig. 5). Phase I is defined as the early response (fast growth) and is followed by a moderate growth phase (Phase II), when the first cohort reaches mesozooplankton compartment. Between $t \sim 30$ and $t \sim 90$ days, the first cohort gradually fades away before a second cohort starts, resulting in a mean mesozooplankton growth close to zero between $t \sim 50$ and $t \sim 80$ days. The time lapse between the two first cohorts defines Phase III. We characterize this phase by three indices: its duration DP3, the day of its onset TP3 and its amplitude AP3 (*ie* the change in growth during the phase). Further details about the calculation of those
<figure><img src="image_12.png" /><figcaption>**Fig. 5** Temporal evolution of the mesozooplankton biomass for all model versions. The time axis is showed on a logarithmic scale. The horizontal lines at the top represent the chronology of the growth phases identified for the 4 models. The dots represent the transition</figcaption></figure>

between two phases, on the chronologies and on the curves. The dotted lines at the end of the chronologies mean that phase V (equilibrium) is not reached before the end of the *HNM* period. The bottom left panel magnifies the model dynamics between days 14 and 18

parameters are provided in supplementary materials (Section 6). For the standard model, DP3 equals 62 days, TP3 equals 61 days and AP3 equals 0.03µmol P L⁻¹d⁻¹. Phase III is then followed by a second moderate growth phase (Phase IV) during which the growth speed gradually decreases until the mesozooplankton biomass reaches an equilibrium at tₐₗ > 1 year (Phase V) (not shown on Fig. 5, discussed in the last section of the results).

**Comparison between the standard and alternative models**

Figure 3 shows that all the model versions (standard and alternative model structures) seem to tend towards the same equilibrium state in *HNM*. Furthermore, the temporal evolution of nutrient concentrations, phytoplankton and

microzooplankton biomasses is similar across all model versions. With the exception of mesozooplankton, on which we focus hereafter, differences between models are weak suggesting that the introduction of mesozooplankton reproduction and size variations in a simple 0-D biogeochemical model such as ours has no significant influence on either lower trophic levels or nutrient availability.

The temporal evolution of the mesozooplankton biomass is displayed on Fig. 5 for all model versions. Overall, the response is slower in the twenty-size-class models (standard and constant-rate models) and in the two-life-stage model than in the no-life-stage model. Consequently, the implementation of reproduction seems to delay the response of mesozooplankton to a sudden increase in food abundance resulting from a phytoplankton bloom and microzooplankton outburst.

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /><figcaption>**Fig. 6** Time evolution of the mesozooplankton sources terms, sinks terms and mean size. All sink and source terms are normalized by the current mesozooplankton biomass *M*. Refer to Eq. 2 for the definition of each term</figcaption></figure>

Using the phases defined for the standard model in the first section of the results, we summarize the differences in the temporal evolution of mesozooplankton between the different model versions (Fig. 5). First, growth is faster in Phase I (initial response) for the constant-rate than for the two-life-stage

and standard models. In the no-life stage model, there is no Phase I. Second, Phase III is not simulated in the no-life-stage and two-life-stage models, so that Phases II and IV are combined into a single pre-equilibrium phase (moderate growth phase). In the constant-rate model, Phase III
cannot be visually identified. Using the indices defined above, we identify however that its temporality is close to the standard model (TP3 = 55 days, DP3 = 55 days) but its amplitude is much lower (AP3 < 0.001 µmolP L⁻¹d⁻¹). Third, Phase V (equilibrium) is not reached before the end of *HNM* for the standard, constant-rate and two-life-stage models (*t*<sub>eq</sub> > 1 year). In the no-life-stage model, Phase V is reached much faster (*t*<sub>eq</sub> ~ 100 days).

Drivers of mesozooplankton dynamics

To explain the differences between the model configurations,
we investigate the different terms of the mesozooplankton
equation (Fig. 6, Eq. 2) and demonstrate the role of explicit
reproduction and of explicit variations in the mean size of
mesozooplankton.

First, we detail the processes involved in the first 15
days (Phase I) in *HNM*. For all models including repro-
duction (two-life-stage, constant rate, standard), the input
from maturation increases sharply to reach 0.5 d⁻¹ during
the second week, and then rapidly decreases to zero at
t = 15 days. This rapid decrease marks the shift from Phase
I to Phase II and happens just after the drop of micro-
zooplankton biomass, which is explained by the sharp
decline in food availability after all phytoplankton have
been grazed by zooplankton. In the size-resolved mod-
els (standard and constant-rate), during the first week of
*HNM*, net growth dominates mesozooplankton dynamics.
Then Phase I actually starts at t = 7 days, and is mostly
explained by the maturation input from the first wave of
juveniles. There is no Phase I in the no-life stage model
because it does not include reproduction.

Second, we detail the processes involved after *t* = 15 days in *HNM* and before the equilibrium is reached. The main difference between the different models is the occurrence of a Phase III in the size-resolved models which is completely absent in the two-life-stage model, and of very low amplitude in the constant-rate model. Phase III is explained by a cohort dynamics characterized by the successive propagation of waves of individuals along the size spectrum. Phase III is actually the period between the first two of these cohorts.

In the standard model, Phase III occurs between *t* ~ 30
days and *t* ~ 90 days and its onset (TP3) corresponds to a
net growth rate minimum. Large size variations in the mean
size of mature organisms during Phase III are shown on the
bottom panel of Fig. 6 but also in Fig. 4. A higher mean size
induces a lower mean maximum growth rate in the standard
model as a consequence of the allometric scaling of growth
rate. Thus, the cohorts materialize in the standard model
dynamics by a strong variation of the net growth rate related
to their ontogenetic growth. In the constant-rate-model, a
slight reduction of the maturation term coincides with Phase
III which occurs slightly earlier than in the standard model

(TP3 = 55 days vs 61 days in the standard model), and has
a slightly lower duration (DP3 = 55 days vs 62 days in the
standard model). A low departure from zero of the growth
closure term coincides with the end of phase III. As growth
rate does not vary in the constant-rate model, net growth
rate variations are lower than in the standard model. Hence,
maturation and growth closure terms dominate the constant-
rate model dynamics during phase III, but the values of the
dominant rates are much lower than in the standard model,
which explains its lower amplitude. Thus, in the constant-
rate model, the cohorts materialize through subtle variations
in maturation and closure terms.

Finally, we analyze the timing of equilibrium (Phase V) for each model version. Mesozooplankton equilibrium is reached once all terms of Eq. 2 are at equilibrium. Differences in this timing are related to the fact that the term reaching equilibrium last differs across models. In the no-life-stage model, Phase V is reached once mortality does not evolve anymore (at *t* ~ 100 days). However, the key process that controls the temporal dynamics of this model is net growth rate, which reaches equilibrium very quickly. Mortality is just following the temporal evolution of the mesozooplankton biomass as it is a density-dependent closure term. In the two-life-stage and constant-rate models, net growth is reaching equilibrium as quickly as in the no-life-stage model. However, maturation continues to evolve until the very end of *HNM* which explains their very long equilibrium time.

As a conclusion on the global behavior of the 4 model
configurations, slower temporal dynamics for all life-stage
models (two-life-stage, constant-rate and standard) than for
the no-life-stage model is induced by the introduction of
reproduction, which converts part of the energy assimilated
by mesozooplankton into juveniles. This results in lower net
growth offset by input from maturation, which introduces a
delay. The appearance of cohort dynamics in models with
explicit size representation is more marked in the case of the
standard model as a consequence of the allometric scaling of
growth rate and explains the faster dynamics in this model
during the first two months.

Sensitivity to the parameters

The biological parameters of our model configurations
have been set to constants and chosen from previous stud-
ies (Table 3). Note that these parameter values are not well
constrained by measurements and observations, either in situ
or in the lab, and that they show large variations within and
between species and groups. The modelling results presented
above may depend on the choice of these parameter val-
ues. To explore the sensitivity of our results to the choice
of parameter values, we performed a sensitivity analysis
<figure><img src="image_15.png" /><figcaption>**Fig. 7** Sensitivity of the mesozooplankton dynamics to the four parameters that have been varied by ±50%. Four indicators of mesozooplankton dynamics are represented: EQ the biomass at equilib-</figcaption></figure>

rium, Q90 the number of days to attain 90% of the mesozooplankton biomass at equilibrium, DP3 the duration of the plateau (in phase III) and TP3 the day of its onset

by systematically varying selected parameters by ±50%. The dynamics of juveniles and mature metazoan organisms show the greatest differences between our four model configuration simulations. Therefore, we focus here on key zooplankton parameters: gross growth efficiency ($\beta$), maximum ingestion rates ($g_Z$ and $g_M$), quadratic mortality rate of mature organisms ($m$) and the fraction of food allocated to reproduction ($w$). Figure 7 shows the changes induced by varying parameter values on five key metrics (mesozooplankton biomass at equilibrium and time to reach equilibrium, and the three metrics used to define Phase III, DP3, TP3, AP3).

The mesozooplankton biomass at equilibrium (EQ) in the 4 model configurations is sensitive to all parameters except that of reproduction. For all parameter values, EQ is similar across the three alternative model versions (constantly rate, two-life-stage and no-life-stage). In the standard model, EQ is slightly larger than in the other model versions for almost

all parameter values. In the supplementary materials (Section 5), we demonstrate this result analytically by manipulating the set of equations of the different model setups.

The time to reach equilibrium (Q90, defined as the number of days to reach 90% of EQ) is also sensitive to all parameters except quadratic mortality. An increase in the maximum ingestion rate or in the growth efficiency leads to a decrease in Q90 for all model configurations because both result in an increase in the mesozooplankton growth rate. When the lowest value of these two parameters is prescribed, mesozooplankton go extinct in the standard, the constant-rate and in the two-life-stage model configurations as the achieved net growth rate remains lower than the dilution rate. But despite the sensitivity to parameter values, the sensitivity analysis confirms one of our main results: the no-life-stage model always simulates the fastest adjustment to equilibrium, while the standard model shows longer Q90 than any other model version, except when using the highest reproduction parameter value.
The cohort dynamics, i.e., the occurrence of Phase III, is simulated in the constant-rate and the standard model configurations for all parameter values, except for the lowest values of the maximum ingestion rate, growth efficiency, and quadratic mortality (for the constant-rate model). In these three cases, a Phase III can not be properly defined, either because mesozooplankton go extinct (the constant-rate model) or because its biomass remains too small (the standard model) to properly identify cohorts. The amplitude AP3 of the Phase III is at least 5 times lower in the constant-rate model than in the standard model for all parameter values. In the standard model, AP3 significantly increases as growth efficiency and microzooplankton maximum growth rate increase. This can be related to the fact that we are considering absolute amplitude, thus directly related to the speed of population growth, which increases as those two parameters increase. The difference in the time of onset TP3 and in the duration DP3 between the constant-rate and the standard model does not exceed 25%, which suggests that the first two cohorts present a similar temporality in both models. DP3 and TP3 systematically decreases when growth efficiency and microzooplankton maximum growth rate increases. Indeed, higher net growth accelerates the dynamics of the cohorts, which are therefore closer together (DP3 decreases) and shorter (TP3 decreases). TP3 and DP3 in the standard model are relatively insensitive to reproduction and quadratic mortality.

# Discussion

## Limitations and prospects for improvement

Our study is based on a series of different model configurations that are used to explore the role of reproduction and ontogenetic growth on the mean state and the temporal dynamics of a plankton system. We employ a very simple 0-D physical framework similar to a chemostat in which the nutrient content of the inflow water is periodically switched between high and low concentrations. This framework has been used in many previous studies either to study the equilibrium states of a modeled system (eg., Thingstad et al. 1996; Ward et al. 2011; Prowe et al. 2012; Grigoratou et al. 2019) or its temporal dynamics, (e.g., Eisenhauer et al. 2009; Banas 2011; Kloosterman et al. 2016). Obviously, chemostat analogues are highly idealized configurations that cannot reproduce the complex and very diverse environmental conditions of the ocean. Furthermore, results are sensitive to the prescribed dilution rate and to the concentrations of the different organisms in the inflow waters that, in particular, controls the viability of the different modeled species and their response to imposed changes (Sommer 1986; Chatterjee and Pal 2013). However, owing to its simplicity, this physical framework represents a useful tool to analyze

and understand the complex behavior of nonlinear systems developed to represent plankton ecosystems.

Our size-resolved description of metazoan organisms does not take into account interspecific and intraspecific differences. The actual body size is the sole structuring trait. Consequently, organisms of the same size are assumed to have the same metabolic and physiological characteristics as well as the same feeding behaviors. Different approaches have been proposed in the literature to introduce diversity into size-spectrum models. Some studies have adopted a mixed formalism between group-or species-based models and size-spectrum models (Blanchard et al. 2009; Maury 2010; Scott et al. 2014). The community is divided into several groups, each being modeled using a size-spectrum approach. A weakness of this mixed formalism is that within each community, body size remains the sole trait and therefore, inter- and intraspecific differences between organisms belonging to the same group remain disregarded. As an alternative to that approach, recent developments have focused on the introduction of a trait-based description of diversity. The most common additional traits are maximum size (Maury and Poggiale 2013) and size at maturity (Hartvig et al. 2011; Castellani et al. 2013). Intraspecific variability is generally represented by the addition of a diffusion term in the McKendrick von Foerster equation (Benoit and Rochet 2004; Datta et al. 2010).

The introduction of a size-resolved representation of mesozooplankton population produces fluctuations of the total biomass of metazoans which correspond to the propagation of successive cohorts through the whole size spectrum (Fig. 4). These fluctuations have been evidenced and extensively studied previously, often using PSPMs, (e.g., McCauley and Murdoch 1987; Persson et al. 1998). These studies showed that these cycles occur when juveniles have a competitive advantage over adults, for instance as a consequence of a higher mass specific ingestion rate (Persson et al. 1998; De Roos and Persson 2003; De Roos et al. 2008; Persson and De Roos 2013). Our size-resolved configurations do assume such an advantage for juveniles as the size scaling factor of ingestion rate is negative. Furthermore, when an allometric relationship is applied to each size class, cyclic cohort dynamics is stronger than when a uniform value is prescribed for juveniles and mature organisms. At equilibrium, cycles disappear despite the competitive advantage of juveniles which seems in contradiction with the previously quoted studies. Our results are difficult to compare to these previous studies which explore the properties of size-structured population models as there are significant differences. For instance, the ressource on which adults and juveniles are able to feed is very often assumed to be the same in many theoretical studies which is not the case in our model experiment: mature organisms feed both on microzooplankton and phytoplankton whereas juveniles only feed on phytoplankton. Furthermore, since our model represents a community,

<figure><img src="image_16.png" /></figure>
organisms can reproduce over a large size range defined by
our mesozooplankton compartment. Finally, a quadratic
mortality rate is applied which differs between the juvenile
and the mature classes. Previous studies have shown that the
way mortality is represented may significantly impact the
behavior of the model, in particular its stability (Murdoch
et al. 2003; van Kooten et al. 2007).

The temporal response of the community size-spectrum
to short time-scales or seasonal perturbations has been inves-
tigated with single size-spectrum approaches that are simi-
lar to ours. They show qualitatively similar results, i.e., the
propagation of successive waves of biomass from small to
large organisms (Pope et al. 1994; Maury et al. 2007; Zhou
et al. 2010). The inclusion of some representation of both
intra- and interspecific diversity has been shown to increase
the intrinsic stability of size-spectrum models (Datta et al.
2010; Zhang et al. 2013). The only study to our knowledge
that has investigated the impact of seasonality with a size-
spectrum model that includes a description of interspecific
biodiversity is the one of Datta and Blanchard (2016). They
showed that the simulated annual mean community size
spectra are identical in seasonal and non seasonal setups. In
addition, each modeled species exhibits propagating waves
from smaller to larger individuals, similarly to what is simu-
lated for the community by single size-spectrum approaches.
Unfortunately, they did not explore the impacts of biodiver-
sity on the community’s response to seasonal fluctuations.
Thus, the impact of both intra- and interspecific diversity on
the response of an ecosystem to changes in the environmen-
tal conditions remains to be studied, at least at intra-seasonal
to seasonal time scales. Therefore, we do not know how the
representation of diversity would alter our results. This could
be investigated in a future study.

The seasonal variability of the environment has resulted
in adaptations in marine organisms, including zooplankton,
in particular at high latitudes, such that their feeding activity,
reproduction and maintenance may display strong seasonal
rhythms (McNamara and Houston 2008). Concerning repro-
duction, our model formulation assumes that zooplankton
are income breeders and thus, allocate part of the available
resources directly to reproduction. This choice was dictated
by simplicity. However, organisms may adopt an alterna-
tive reproduction strategy which is termed capital breed-
ing. During the feeding season, an individual may allocate
energy to reserves that are used later in the year to reproduce
(Jönsson 1997; Varpe et al. 2009). Capital breeding has been
suggested to be a better strategy in highly seasonal environ-
ments with a short feeding season, such as the polar oceans
(Sainmont et al. 2014). Obviously, changing the breeding
mode in our model configuration would significantly change
the simulated temporal dynamics of juveniles and mature
organisms as reproduction would become temporally decou-
pled from feeding.

Implications

Mesozooplankton occupy a key position in most ocean food chains. As predators of smaller plankton, they can influence the dynamics of microzooplankton and phytoplankton. As producers of fecal pellets (Steinberg and Landry 2017), and since they are able to migrate vertically in the water column (Emerson et al. 1997), they also control the transfer of matter from the ocean surface to depth, and as such may influence the marine carbon and nutrient cycles. As preys for larger organisms, mesozooplankton also control the transfer of energy to the higher trophic levels (Beaugrand et al. 2010). By eating and being eaten, mesozooplankton both exert a top-down control (on marine biogeochemistry and lower-trophic levels) and a bottom-up control (on higher trophic levels) (Kiørboe 1997). Bottom-up effects induced by the size structuration of a population with size-dependent life-history traits have been summarized by De Roos (2020). In particular, a predator feeding on a size structured prey may suffer from an Allee effect (i.e., reduced reproductive success at low population densities) with the emergence of a stable state where the predator population collapses, a behavior that is not observed in the absence of size structuring (De Roos and Persson 2002).

Here, we show that the explicit representation of mea-
tozoans reproduction and size-growth in a idealized 0-D
biogeochemical model strongly affect the temporal dynam-
ics of mesozooplankton. However, we do not highlight any
effects on the temporal dynamics of phytoplankton and
nutrients when including reproduction and growth in size.
Previous studies have on the contrary demonstrated
potential top-down effects. It has been shown that
early life stages of copepods, such as nauplii, can represent a
significant part of microzooplankton (up to 30%) (Quevedo
and Anadón 2000; Safi et al. 2007) and hence induce a top
down control of mesozooplankton abundance (in this case
with the species M. leidyi) on microzooplankton (McNamara
et al. 2013). Also, mesozooplankton have been shown to be
responsible for a significant part of the vertical export of
carbon due to fecal pellets (Stamieszkin et al. 2015). Explor-
ing those latter processes implies a vertically resolved rep-
resentation of biogeochemical processes, which is absent in
our model. This could explain why we do not observe any
effect of our different representations of mesozooplankton on
nutrients and lower-trophic levels.

As an explicit representation of individual size growth
and reproduction of mesozooplankton in our model induces
a delay in the mesozooplankton response to a phytoplankton
bloom, we could expect that in a model coupled with upper
trophic levels, the response of mesozooplankton predators
would also be delayed. Indeed, numerous examples of bot-
tom-up effect of mesozooplankton on upper trophic levels
have been described which calls to the developments of
robust zooplankton submodels in end-to-end marine ecosystem models (Rose et al. 2010; Mitra et al. 2014). As copepods dominate mesozooplankton in many systems (Gallienne 2001), studies on those species may give a representative overview on the bottom-up control of mesozooplankton on upper trophic levels. Bi et al. (2011) focused on the correlation between water-type affinities of copepod communities and salmon type predominancy along the Washington and Oregon coast and have showed that the predominant salmon species shifts followed switchs between cold- and warm-water copepod species, which suggest that salmon population dynamics is controlled by prey availability. In the Mediterranean Sea, Saraux et al. (2019) have shown that variations in plankton community around 2008, especially marked by a decrease in copepod abundance, have had a direct impact on fish biomass and on the mean size and weight of sardine and anchovy, highlighting again an important bottom-up effect.

## Conclusion

The explicit representation of metazoan life-cycle in our simple biogeochemical model does not affect the temporal dynamics of non-mesozooplankton compartments. However, mesozooplankton temporal dynamics is changed, mostly through a few months delay in the response to an increase in food availability when reproduction is taken into account. In model versions with reproduction, we also identify an early acceleration phase, in which rates are particularly high, but which happens before mesozooplankton biomass reaches a significant value. Lastly, accounting for ontogenetic growth induces cohort dynamics in metazoans. The comparison of our 4 model versions enables to get some insights into the main factors inducing these changes.

First, as a consequence of size structuration, the temporal dynamics is affected by the size-based parameterization of the maximal growth rates. At the start of HNM, when zooplankton biomass is limited, metazoan individuals are mostly juveniles, i.e., they are part of the microzooplankton and therefore benefit from high growth rates. Those fast growing individuals are then matured into adult animals, inducing a faster initial increase of mesozooplankton to nutrient input as compared to a model without reproduction. In addition, temporal changes in the size distribution of mesozooplankton lead to the appearance of cohort dynamics, totally absent in model configurations without discretization of metazoan size.

The second factor influencing the temporal dynamics is the introduction of non-predation terms between microzooplankton and mesozooplankton size-classes, i.e., reproduction and maturation. Differences in the temporal dynamics of those two opposite terms induce a delay in the mesozooplankton response to equilibrium: reproduction, a loss term

for mature mesozooplankton, reaches its equilibrium faster than the gain term, maturation. This results in a net loss term and reduces mesozooplankton growth before equilibrium is reached, and hence retards this equilibrium when compared to a model without any explicit representation of reproduction.

Mesozooplankton play a key role in marine ecosystems because of their intermediate position in trophic food webs. Explicit representation of mesozooplankton reproduction and ontogenetic growth should be considered when modeling the upper trophic levels as the dynamics of upper trophic levels feeding on mesozooplankton could be affected, particularly in regions where the seasonality is marked. Numerous top down effect on lower trophic levels (particularly microzooplankton) and biogeochemistry (particularly on the vertical carbon export), not testable in the context of our 0-D study, would be expected in a vertically resolved model. Size spectrum approaches are interesting avenues that could provide all the tools necessary to improve mesozooplankton community representation in biogeochemical models.

**Supplementary information** The online version of this article (https://doi.org/10.1007/s12080-021-00519-5) contains supplementary material, which is available to authorized users.

**Acknowledgements** We thank Quentin Jutard and Ben Mayersohn for providing a basis framework to develop our model. We are grateful to the two anonymous reviewers who have helped us improve the manuscript.

**Author contributions** CC, OA and LB contributed to the study conception and design. Model development and analysis were performed by CC. All authors contributed to the interpretation of the results. CC wrote a first draft of the manuscript, which was then edited by all authors. All authors read and approved the final manuscript.

**Funding** This study has received funding from the Agence Nationale de la Recherche grant agreement ANR-17-CE32-0008 (CIGOEF).

**Availability of data and materials** Codes for the model presented in this paper are available on Zenodo with the identifier https://doi.org/10.5281/zenodo.5040857.

## Declarations

**Consent for publication** All authors read and approved the final manuscript.

**Competing interests** The authors declare that they have no conflict of interest.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in
the article's Creative Commons licence and your intended use is not
permitted by statutory regulation or exceeds the permitted use, you will
need to obtain permission directly from the copyright holder. To view a
copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

References

Aumont O, Ethé C, Tagliabue A, Bopp L, Gehlen M (2015) PISCES-v2: An ocean biogeochemical model for carbon and ecosystem studies. Geosci Model Dev 8(8):2465-2513

Aumont O, Maury O, Lefort S, Bopp L (2018) Evaluating the Potential Impacts of the Diurnal Vertical Migration by Marine Organisms on Marine Biogeochemistry. Global Biogeochem Cycles 32(11):1622-1643

Baird ME, Suthers IM (2007) A size-resolved pelagic ecosystem model. Ecol Model 203(3-4):185-203

Banas NS (2011) Adding complex trophic interactions to a size-spectral plankton model: Emergent diversity patterns and limits on predictability. Ecol Model 222(15):2663-2675

Basedow SL, Tande KS (2006) Cannibalism by female Calanus finmarchicus on naupliar stages. Mar Ecol Prog Ser 327:247-255

Beaugrand G, Edwards M, Legendre L (2010) Marine biodiversity, ecosystem functioning, and carbon cycles. Proc Natl Acad Sci 107(22):10120-10124

Benoit E, Rochet MJ (2004) A continuous model of biomass size spec-
tra governed by predation and the effects of fishing on them. J
Theor Biol 226(1):9-21

Bi H, Peterson WT, Lamb J, Casillas E (2011) Copepods and salmon: Characterizing the spatial distribution of juvenile salmon along the Washington and Oregon coast, USA. Fish Oceanogr 20(2):125-138

Bissinger JE, Montagnes DJ, Sharples J, Atkinson D (2008) Predicting marine phytoplankton maximum growth rates from temperature: Improving on the Eppley curve using quantile regression. Limnol Oceanogr 53(2):487-493

Blanchard JL, Jennings S, Law R, Castle MD, McCloghrie P, Rochet M-J, Benot E (2009) How does abundance scale with body size in coupled size-structured food webs? J Anim Ecol 78(1):270-280

Brodeur RD, Terazaki M (1999) Springtime abundance of chaetognaths in the shelf region of the northern Gulf of Alaska, with observations on the vertical distribution and feeding of Sagitta elegans. Fish Oceanogr 8(2):93-103

Buitenhuis E, Qur CL, Aumont O, Beaugrand G, Bunker A, Hirst A, Ikeda T, O'Brien T, Piontkovski S, Straile D (2006) Biogeochemical fluxes through mesozooplankton. Global Biogeochem Cycles 20:2

Calbet A (2001) Mesozooplankton Grazing Effect on Primary Pro-
duction: A Global Comparative Analysis in Marine Ecosystems.
Limnol Oceanogr 46:1824-1830

Calbet A, Landry MR (2004) Phytoplankton growth, microzooplankton grazing, and carbon cycling in marine systems. Limnol Oceanogr 49(1):51-57

Carlotti F, Giske J, Werner F (2000) 12 - Modeling zooplankton dynamics. In: Harris R, Wiebe P, Lenz J, Skjoldal HR, Huntley M (eds) ICES Zooplankton Methodology Manual. Academic Press, London, pp 571-667

Carlotti F, Poggiale JC (2010) Towards methodological approaches to implement the zooplankton component in end to end food-web models. Prog Oceanogr 84(1):20-38

Castellani M, Rosland R, Urtizberea A, Fiksen Ø (2013) A mass-balanced pelagic ecosystem model with size-structured behaviourally adaptive zooplankton and fish. Ecol Model 251:54-63

Chatterjee A, Pal S (2013) Coexistence of plankton model with essential multiple nutrient in chemostat. Int J Biomath 06(06):1350042

Chisholm SW (1992) Phytoplankton Size. Springer, US, Boston, MA, pp 213–237

Datta S, Blanchard JL (2016) The effects of seasonal processes on size spectrum dynamics. Can J Fish Aquat Sci 73(4):598–610

Datta S, Delius GW, Law R (2010) A Jump-Growth Model for Predator Prey Dynamics: Derivation and Application to Marine Ecosystems. Bull Math Biol 72(6):1361–1382

De Roos AM (1997) A gentle introduction to physiologically struc-
tured population models. In Structured-population models
in marine, terrestrial, and freshwater systems. Springer, pp.
119–204

De Roos AM (2018) When individual life history matters: conditions for juvenile-adult stage structure effects on population dynamics. Thyroid Res 11(4):397-416

De Roos AM (2020) Effects of life history and individual development on community dynamics: A review of counterintuitive consequences. Ecol Res 35(6):930-946

De Roos AM, Persson L (2002) Size-dependent life-history traits promote catastrophic collapses of top predators. Proc Natl Acad Sci 99(20):12907-12912

De Roos AM, Persson L (2003) Competition in size-structured popu-
lations: Mechanisms inducing cohort formation and population
cycles. Theor Popul Biol 63(1):1-16

De Roos AM, Persson L, McCauley E (2003) The influence of size-
dependent life-history traits on the structure and dynamics of
populations and communities. Ecol Lett 6(5):473-487

De Roos AM, Schellekens T, Van Kooten T, Van De Wolfshaar K, Claessen D, Persson L (2008) Simplifying a physiologically structured population model to a stage-structured biomass model. Theor Popul Biol 73(1):47-62

Edwards KF, Thomas MK, Klausmeier CA, Litchman E (2012) Allo-
metric scaling and taxonomic variation in nutrient utilization
traits and maximum growth rate of phytoplankton. Limnol.
Oceanogr 57(2):554-566

Edwards KF, Thomas MK, Klausmeier CA, Litchman E (2016) Phyto-
plankton growth and the interaction of light and temperature: A
synthesis at the species and community level. Limnol Oceanogr
61(4):1232-1244

Eisenhauer L, Carlotti F, Baklouti M, Diaz F (2009) Zooplankton population model coupled to a biogeochemical model of the North Western Mediterranean Sea ecosystem. Ecol Model 220(21):2865-2876

Emerson S, Quay P, Karl D, Winn C, Tupas L, Landry M (1997) Experimental determination of the organic carbon flux from open-ocean surface waters. Nature 389(6654):951-954

Eppley RW (1972) Temperature and phytoplankton growth in the sea. Fish. bull 70(4):1063-1085

Eppley RW, Rogers JN, McCarthy JJ (1969) Half-saturation constants for uptake of nitrate and ammonium by marine phytoplankton. Limnol Oceanogr 14(6):912-920

Evans GT, Parslow JS, Evans GT, Parslow JS (1985) A model of annual plankton cycles. Deep Sea Research Part B. Oceanogr Lit Rev 32(9):759

Everett JD, Baird ME, Buchanan P, Bulman C, Davies C, Downie R, Griffiths C, Heneghan R, Kloser RJ, Laiolo L, Lara-Lopez A, Lozano-Montes H, Matear RJ, McEnnulty F, Robson B, Rochester W, Skerratt J, Smith JA, Strzelecki J, Suthers IM, Swadling KM, van Ruth P, Richardson AJ (2017) Modeling What We Sample and Sampling What We Model: Challenges for Zooplankton Model Assessment. Front Mar Sci 4
Fasham MJR, Ducklow HW, McKelvie SM (1990) A nitrogen-based model of plankton dynamics in the oceanic mixed layer. J Mar Res 48(3):591–639

Fennel W (2001) Modeling of copepods with links to circulation models. J Plankton Res 23(11):1217–1232

Flynn KJ, Irigoien X (2009) Aldehyde-induced insidious effects cannot be considered as a diatom defence mechanism against copepods. Mar Ecol Prog Ser 377:79–89

Follows MJ, Dutkiewicz S, Grant S, Chisholm SW (2007) Emergent biogeography of microbial communities in a model ocean. Science 315(5820):1843–1846

Gallienne CP (2001) Is Oithona the most important copepod in the world's oceans? J Plankton Res 23(12):1421–1432

Gentleman W, Leising A, Frost B, Strom S, Murray J (2003) Functional responses for zooplankton feeding on multiple resources: a review of assumptions and biological dynamics. Deep-Sea Res II(50):2847–2875

Gentleman WC, Neuheimer AB (2008) Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, food-limitation and acclimation. J Plankton Res 30(11):1215–1231

Grigoratou M, Monteiro FM, Schmidt DN, Wilson JD, Ward BA, Ridgwell A (2019) A trait-based modelling approach to planktonic foraminifera ecology. Biogeosciences 16:1469–1492

Guiet J, Poggiale JC, Maury O (2016) Modelling the community size-spectrum: recent developments and new directions. Ecol Model 337:4–14

Guill C (2009) Alternative dynamical states in stage-structured consumer populations. Theor Popul Biol 76(3):168–178

Hansen PJ, Bjørnsen PK, Hansen BW (1997) Zooplankton grazing and growth: Scaling within the 2–2,000-µm body size range. Limnol Oceanogr 42(4):687–704

Hartvig M, Andersen KH, Beyer JE (2011) Food web framework for size-structured populations. J Theor Biol 272(1):113–122

Henschke N, Everett JD, Richardson AJ, Suthers IM (2016) Rethinking the role of salps in the ocean. Trends Ecol Evol 31(9):720–733

Hernandez-Len S, Ikeda T (2005) A global assessment of mesozooplankton respiration in the ocean. J Plankton Res 27(2):153–158

Jönsson KI (1997) Capital and Income Breeding as Alternative Tactics of Resource Use in Reproduction. Oikos 78(1):57–66

Kiørboe T (1997) Population regulation and role of mesozooplankton in shaping marine pelagic food webs. Hydrobiologia 363(1):13–27

Kloosterman M, Campbell SA, Poulin FJ (2016) An NPZ Model with State-Dependent Delay Due to Size-Structure in Juvenile Zooplankton. SIAM J Appl Math 76(2):551–577

Kooijman SALM (2013) Waste to hurry: dynamic energy budgets explain the need of wasting to fully exploit blooming resources. Oikos 122(3):348–357

Le Queré C, Harrison, SP, Colin PI, Buitenhuis ET, Aumont O, Bopp L, Claustre H, Cotrim Da Cunha L, Geider R, Giraud X, Klaas C, Kohfeld KE, Legendre L, Manizza M, Platt T, Rivkin RB, Sathyendranath S, Uitz J, Watson AJ, Wolf-Gladrow D (2005) Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models. Glob Chang Biol

Maran E, Cermeo P, Lopez-Sandoval DC, Rodriguez-Ramos T, Sobrino C, Huete-Ortega M, Blanco JM, Rodriguez J (2013) Unimodal size scaling of phytoplankton growth and the size dependence of nutrient uptake and use. Ecol Lett 16(3):371–379

Maury O (2010) An overview of APECOSM, a spatialized mass balanced Apex Predators ECOSystem Model to study physiologically structured tuna population dynamics in their ecosystem. Prog Oceanogr 84(1–2):113–117

Maury O, Poggiale JC (2013) From individuals to populations to communities: A dynamic energy budget model of marine ecosystem size-spectrum including life history diversity. J Theor Biol 324:52–71

Maury O, Shin YJ, Faugeras B, BenAri T, Marsac F (2007) Modeling environmental effects on the size-structured energy flow through marine ecosystems. Part 2: Simulations. Prog Oceanogr 74(4):500–514

Mayersohn B (2018) Phyre. https://github.com/benmayersohn/phyre

McCauley E, Murdoch WW (1987) Cyclic and stable populations: Plankton as paradigm. Am Nat 129(1):97–121

McNamara JM, Houston AI (2008) Optimal annual routines: behaviour in the context of physiology and ecology. Philosophical Transactions of the Royal Society B: Biological Sciences 363(1490):301–319

McNamara ME, Lonsdale DJ, Cerrato RM (2013) Top-down control of mesozooplankton by adult mnemiopsis leidyi influences microplankton abundance and composition enhancing prey conditions for larval ctenophores. Estuar Coast Shelf Sci 133:2–10

Metz J, Diekmann O (1986) A gentle introduction to structured population models: three worked examples. In The dynamics of physiologically structured populations. Springer, pp.3–45

Mitra A, Castellani C, Gentleman WC, Jnasdtir SH, Flynn KJ, Bode A, Halsband C, Kuhn P, Licandro P, Agersted MD, Calbet A, Lindeque PK, Koppelmann R, Møller EF, Gislason A, Nielsen TG, St. John M (2014) Bridging the gap between marine biogeochemical and fisheries sciences; configuring the zooplankton link. Prog Oceanogr 129, Part B, 176–199

Moore JK, Doney SC, Kleypas JA, Glover DM, Fung IY (2002) An intermediate complexity marine ecosystem model for the global domain. Deep Sea Res II(49):403–462

Murdoch WW, Briggs CJ, Nisbet RM (2003) Consumer-resource dynamics, vol. 36. Princeton University Press

Ohman MD, Durbin EG, Runge JA, Sullivan BK, Field DB (2008) Relationship of predation potential to mortality of Calanus finmarchicus on Georges Bank, northwest Atlantic. Limnol Oceanogr 53(4):1643–1655

Packard TT, Gomez M (2013) Modeling vertical carbon flux from zooplankton respiration. Prog Oceanogr 110:59–68

Persson L, De Roos AM (2013) Symmetry breaking in ecological systems through different energy efficiencies of juveniles and adults. Ecology 94(7):1487–1498

Persson L, Leonardsson K, De Roos AM, Gyllenberg M, Christensen B (1998) Ontogenetic scaling of foraging rates and the dynamics of a size-structured consumer-resource model. Theor Popul Biol 54(3):270–293

Pinceel T, Vanschoenwinkel B, Brendonck L, Buschke F (2016) Modelling the sensitivity of life history traits to climate change in a temporary pool crustacean. Sci Rep 6(1):1–5

Pope JG, Shepherd JG, Webb J, Stebbing ARD, Mangel M, Beverton RJH, Beddington JR, Beddington JR, Cushing DH, May RM, Steele JH (1994) Successful surf-riding on size spectra: the secret of survival in the sea. Philosophical Transactions of the Royal Society of London. Series B: Biological Sciences 343(1303):41–49

Prowe AEF, Pahlow M, Oschlies A (2012) Controls on the diversity productivity relationship in a marine ecosystem model. Ecol Model 225:167–176

Quevedo M, Anadón R (2000) Spring microzooplankton composition, biomass and potential grazing in the central Cantabrian coast (southern Bay of Biscay). Oceanol Acta 23(3):297–310

Quinones RA, Platt T, Rodriguez J (2003) Patterns of biomass-size spectra from oligotrophic waters of the Northwest Atlantic. Prog Oceanogr 57(3):405–427

Rose KA, Allen JI, Artioli Y, Barange M, Blackford J, Carlotti F, Cropp R, Daewel U, Edwards K, Flynn K, Hill SL, HilleRisLambers R, Huse G, Mackinson S, Megrey B, Moll A, Rivkin R, Salihoglu B, Schrum C, Shannon L, Shin YJ, Smith SL, Smith C, Solidoro C, St.John M, Zhou M (2010) End-to-end models for the analysis of marine ecosystems: Challenges, issues, and next steps. Marine and Coastal Fisheries 2(1):115–130
Safi KA, Brian Griffiths F, Hall JA (2007) Microzooplankton com-
position, biomass and grazing rates along the WOCE SR3 line
between Tasmania and Antarctica. Deep-Sea Res I Oceanogr Res
Pap 54(7):1025-1041

Sailley S, Vogt M, Doney S, Aita M, Bopp L, Buitenhuis E, Hashioka
T, Lima I, Qur CL, Yamanaka Y (2013) Comparing food web
structures and dynamics across a suite of global marine ecosys-
tem models. Ecol Model 261-262:43-57

Sainmont J, Andersen KH, Varpe Ø, Visser AW (2014) Capital versus Income Breeding in a Seasonal Environment. Am Nat 184(4):466-476

Saraux C, Beveren EV, Brosset P, Queiros Q, Bourdeix JH, Dutto G, Gasset E, Jac C, Bonhommeau S, Fromentin JM (2019) Small pelagic fish dynamics: A review of mechanisms in the gulf of lions. Deep Sea Research Part II: Topical Studies in Oceanography. Drivers of dynamics of small pelagic fish resources: environmental control of long-term changes 159:52-61

Schmoker C, Hernandez-Len S, Calbet A (2013) Microzooplankton grazing in the oceans: impacts, data variability, knowledge gaps and future directions. J Plankton Res 35(4):691-706

Scott F, Blanchard JL, Andersen KH (2014) mizer: an R package for multispecies, trait-based and community size spectrum ecological modelling. Methods Ecol Evol 5(10):1121-1125

Serra-Pompei C, Soudijn F, Visser AW, Kiørboe T, Andersen KH (2020) A general size- and trait-based model of plankton communities. Prog Oceanogr 189, August, 102473

Sheldon RW, Prakash A, Sutcliffe WH (1972) The size distribution of particles in the ocean. Limnol Oceanogr 17(3):327-340

Sommer U (1986) Phytoplankton competition along a gradient of dilution rates. Oecologia 68(4):503-506

Soudijn FH, De Roos AM (2017) Approximation of a physiologically structured population model with seasonal reproduction by a stage-structured biomass model. Thyroid Res 10(1):73-90

Stamieszkin K, Pershing AJ, Record NR, Pilskaln CH, Dam HG, Feinberg LR (2015) Size as the master trait in modeled copepod fecal pellet carbon flux. Limnol Oceanogr 60(6):2090-2107

Steinberg DK, Landry MR (2017) Zooplankton and the ocean carbon cycle. Ann Rev Mar Sci 9(1):413-444 (PMID: 27814033)

Stock C, Dunne J (2010) Controls on the ratio of mesozooplankton production to primary production in marine ecosystems. Deep-Sea Res I Oceanogr Res Pap 57(1):95-112

Straile D (1997) Gross growth efficiencies of protozoan and metazoan zooplankton and their dependence on food concentration, predator-prey weight ratio, and taxonomic group. Limnol Oceanogr 42(6):1375-1385

Thingstad TF, Havskum H, Garde K, Riemann B (1996) On the strategy of 'eating your competitor': A mathematical analysis of algal mixotrophy. Ecology 77(7):2108-2118

Thompson GA, Dinofrio EO, Alder VA (2013) Structure, abundance and biomass size spectra of copepods and other zooplankton communities in upper waters of the Southwestern Atlantic Ocean during summer. J Plankton Res 35(3):610-629

van Kooten T, Persson L, De Roos AM (2007) Size-dependent mortality induces life-history changes mediated through population dynamical feedbacks. Am Nat 170(2):258-270

Varpe Ø, Jørgensen C, Tarling GA, Fiksen Ø (2009) The adaptive value of energy storage and capital breeding in seasonal environments. Oikos 118(3):363-370

Ward BA, Dutkiewicz S, Barton AD, Follows MJ (2011) Biophysical Aspects of Resource Acquisition and Competition in Algal Mixotrophs. Am Nat 178(1):98-112

Ward BA, Dutkiewicz S, Jahn O, Follows MJ (2012) A size-structured food-web model for the global ocean. Limnol Oceanogr 57(6):1877-1891

Wirtz KW (2011) Non-uniform scaling in phytoplankton growth rate due to intracellular light and CO2 decline. J Plankton Res 33(9):1325-1341

Zhang L, Thygesen UH, Knudsen K, Andersen KH (2013) Trait diversity promotes stability of community dynamics. Thyroid Res 6(1):57-69

Zhou M, Carlotti F, Zhu Y (2010) A size-spectrum zooplankton closure model for ecosystem modelling. J Plankton Res 32(8):1147-1165

Zhou M, Tande KS, Zhu Y, Basedow S (2009) Productivity, trophic levels and size spectra of zooplankton in northern Norwegian shelf regions. Deep-Sea Res II Top Stud Oceanogr 56(21):1934-1944