# On the sensitivity of plankton ecosystem models to the formulation of zooplankton grazing

Fanny Chenillat, Pascal Rivière, Mark Ohman

► To cite this version:

Fanny Chenillat, Pascal Rivière, Mark Ohman. On the sensitivity of plankton ecosystem models to the formulation of zooplankton grazing. PLoS ONE, 2021, 16 (5), pp.e0252033. 10.1371/journal.pone.0252033. hal-03324377

HAL Id: hal-03324377

https://hal.univ-brest.fr/hal-03324377

Submitted on 27 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
RESEARCH ARTICLE

# On the sensitivity of plankton ecosystem models to the formulation of zooplankton grazing

Fanny Chenillat<sup>1*</sup>, Pascal Rivière<sup>1</sup>, Mark D. Ohman<sup>2</sup>

1 Laboratoire des Sciences de l'Environnement Marin (LEMAR), Université de Brest, Ifremer, IRD, IUEM, Brest, France, 2 Scripps Institution of Oceanography, University of California, San Diego, La Jolla, California, United States of America

* fanny.chenillat@univ-brest.fr

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

## Abstract

Model representations of plankton structure and dynamics have consequences for a broad spectrum of ocean processes. Here we focus on the representation of zooplankton and their grazing dynamics in such models. It remains unclear whether phytoplankton community composition, growth rates, and spatial patterns in plankton ecosystem models are especially sensitive to the specific means of representing zooplankton grazing. We conduct a series of numerical experiments that explicitly address this question. We focus our study on the form of the functional response to changes in prey density, including the formulation of a grazing refuge. We use a contemporary biogeochemical model based on continuum size-structured organization, including phytoplankton diversity, coupled to a physical model of the California Current System. This region is of particular interest because it exhibits strong spatial gradients. We find that small changes in grazing refuge formulation across a range of plausible functional forms drive fundamental differences in spatial patterns of plankton concentrations, species richness, pathways of grazing fluxes, and underlying seasonal cycles. An explicit grazing refuge, with refuge prey concentration dependent on grazers' body size, using allometric scaling, is likely to provide more coherent plankton ecosystem dynamics compared to classic formulations or size-independent threshold refugia. We recommend that future plankton ecosystem models pay particular attention to the grazing formulation and implement a threshold refuge incorporating size-dependence, and we call for a new suite of experimental grazing studies.

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

### OPEN ACCESS

**Citation:** Chenillat F, Rivière P, Ohman MD (2021) On the sensitivity of plankton ecosystem models to the formulation of zooplankton grazing. PLoS ONE 16(5): e0252033. https://doi.org/10.1371/journal.pone.0252033

**Editor:** Elena Gorokhova, Stockholm University, SWEDEN

**Received:** February 14, 2021

**Accepted:** May 8, 2021

**Published:** May 25, 2021

**Copyright:** © 2021 Chenillat et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All model data files are available from the figshare database (accession number https://doi.org/10.6084/m9 figshare. 12696422.v1). All observation data files are available online: oceaninformatics.ucsd.edu/datazoo.

**Funding:** This work was supported by four grants: (1) F.C. received a LabexMer grant ANR-10-LABX-19-01 (www.labexmer.eu), (2) F.C. received a the PRESTIGE-2015-3-0017 grant from the People Programme (Marie Curie Actions - PRESTIGE program coordinated - by Campus France via the

## Introduction

Marine plankton dynamics are key determinants of ocean processes ranging from biogeochemical cycling [1], sustainable fisheries [2], and persistence of endangered and protected species, to local [3] and global [4,5] patterns of ocean biodiversity. Phytoplankton growth can alter ocean heat budgets [6]), and zooplankton behavior can potentially modify turbulent mixing [7]. Hence, model representations of plankton structure and dynamics have consequences
European Union's Seventh Framework Programme (FP7/2007-2013) under REA grant agreement n. PCOFUND-GA-2013-609102) (www.campusfrance.org); (3) M.D.O. received support from the Gordon and Betty Moore Foundation (Grant id GBMF3576 and GBMF5479) (www.moore.org) (4) M.D.O. received a U.S. NSF grants to the California Current Ecosystem LTER site (#1637632) (www.nsf.gov/) The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

for a broad spectrum of contemporary issues. Yet, representing the structure of plankton eco-
systems in tractable ocean models can be daunting, since holoplanktonic organisms include
several thousand species [8,9], untold numbers of Operational Taxonomic Units, at least six
trophic levels, and widespread mixotrophy [10]. Moreover, such organisms span a size range
(as biovolume) of at least 16 orders of magnitude, a particularly relevant point because preda-
tor-prey interactions are thought to be at least partially size-dependent [11]. Approaches taken
to simplify representation of this taxonomic, trophic, and size diversity range from early NPZ
(nutrient-phytoplankton-zooplankton) models [12–14] that compressed all dynamics to three
internally uniform trophic levels, to plankton functional types [15–19] and trait-based
approaches (Martini et al. 2020 [20] and references within). Some have represented phyto-
plankton size classes with allometrically scaled relationships [21–25], while ‘self-organizing’
models permit emergence of phytoplankton types of different physiological characteristics
[4,26]. Among these approaches, allometric scaling has been appealing to many because the
number of free parameters in a model can be dramatically reduced.

Since many important scientific issues require analysis on large spatial scales, plankton eco-
system models are increasingly embedded in Global Circulation Models. This scale is neces-
sary to address, for example, ocean carbon sequestration partly mediated by the planktonic
food web in the Biological Carbon Pump and sustainable fisheries production. However, to
date such global scale models have been more successful at reproducing primary production
and spatial distributions of plankton types in the open ocean [4,24,25,27] than in the dynamic
boundary currents at the eastern and western margins of the ocean basins. Eastern Boundary
Currents, in particular, are regions of elevated and highly dynamic carbon fluxes [28,29] and
fisheries production [30], suggesting that plankton growth and composition in such regions
need to be much better represented than currently found in global models.

Although it is common to explain the deficiencies of biogeochemical model fits by pointing
to the more complex physics and need for higher spatial resolution in Eastern Boundary Curr-
ents [29,31], here we draw attention, instead, to the critical importance of more accurate
representation of biological processes. We focus, in particular, on the zooplankton and their
grazing dynamics. Historically, much of the development of plankton ecosystem modeling has
emphasized phytoplankton growth, with the assumption that phytoplankton temporal dynam-
ics and spatial patterns are governed primarily by responses to temperature, light, and nutri-
ents (including trace metals). The effects of these variables on phytoplankton growth and
species composition are generally well constrained by laboratory experiments and readily
incorporated into model structures (e.g., [4,18,24,32]), although interaction terms among vari-
ables need further attention (e.g., [33]). In contrast to controls on the specific growth rate of
phytoplankton, loss terms, principally grazing by zooplankton [34] but also viral lysis and sink-
ing, are typically represented in a highly simplified form. In many respects, this imbalance is
surprising, since the functional response of zooplankton grazing to changes in prey density is
known to markedly affect phytoplankton dynamics (e.g., [12,27,35–40]) and top-down grazing
control has been suggested to play a fundamental role in productive ecosystems [24]. In addi-
tion, the choice of mathematical formulation of grazing, including subtle differences, has been
shown to induce large deviations in the plankton model response [12,27,41,42]. Nevertheless,
the imbalance of effort and attention to growth processes over grazing losses persists in the lit-
erature, especially for models applied on the large scale (cf. [4]).

A consequence of this imbalanced effort is that it remains unclear whether the specific
means of representing zooplankton grazing is likely to have a governing role for phytoplank-
ton community composition, growth rates, and spatial patterns in contemporary, spatially-
resolved plankton ecosystem models. Here we conduct a series of numerical experiments that
explicitly address this question in order to assess the sensitivity of phytoplankton spatial
distributions, community composition, and food web energy flow to the formulation of zoo-
plankton grazing. We do not seek to explore detailed taxon-specific grazing behaviors (e.g.,
[43]), which are currently computationally intractable for global-scale models, nor do we
attempt to expand the number and diversity of types of modeled zooplankton represented.
Instead, we focus on a universal aspect of phytoplankton-grazer interactions, i.e., the form of
the functional response to changes in prey density, including the formulation of a grazing ref-
uge [41,42]. We evaluate the consequences that alternative representations may have for key
ecosystem processes. Anderson et al. [27] analyzed the consequences of four different func-
tional responses for plankton spatial distributions, primary production, and carbon export,
using two zooplankton groups, but did not explicitly test sensitivity to grazing refugia.

We conduct a process-oriented numerical study with a size-structured plankton ecosystem model that has been previously calibrated and validated on a global scale [24]. We perform the numerical simulations with a physical model of the California Current System (CCS), an Eastern Boundary Current upwelling ecosystem where better representation of food web processes is needed. In the context of this process study, we analyze the response of the ecosystem in terms of structure and functioning by varying only the form of the grazing term. We are thus able to quantify the changes induced by such process and we show that an explicit grazing refugium, with refuge prey concentration dependent on grazers’ body size, is likely to provide more coherent plankton ecosystem dynamics compared to more classical formulations as in the original model.

**Materials and methods**

In this study we adopt a process study methodology to focus on the effects of grazing formula-
tion on ecosystem structure and functioning. We use a physical regional model that was previ-
ously validated in the California Current System to force a continuum size-structured
ecosystem model. The ecosystem model was previously calibrated and validated on global
scale. Our goal is not to reproduce the real ecosystem dynamics by calibrating or finely tuning
the model to specific field conditions, but rather to study its sensitivity to the grazing formula-
tion. This is why we only change the grazing term of the biological model equations; all other
terms and parameters being kept to their original forms and values. We describe and quantify
the consequences of different grazing formulations on the ecosystem dynamics in terms of
Chl-a distribution, seasonal plankton biomass evolution, trophic interactions and plankton
spatial niches and richness.

**Model configuration**

**Physical model.** We use the Regional Ocean Modeling System (ROMS), a hydrodynamic model that is a three dimensional, free-surface, hydrostatic, eddy-resolving primitive equation ocean model based on the Boussinesq approximation and hydrostatic vertical momentum balance [44]. The model configuration is the same that Chenillat et al. [45] successfully used with an NPZD-type model over the entire CCS, from the coast to 1000 km offshore and from Baja California (24°N) to Vancouver Island (50°N). It has a 15 km horizontal resolution with 32 sigma-coordinate levels irregularly spaced with higher resolution near the surface to adequately resolve upper ocean physics and ecosystem dynamics. We use monthly mean climatologies to force the model at the surface and its lateral boundaries, as in Chenillat et al. [46,47].

After a 10-year physical spin-up, the coupled biological-physical model is integrated for-
ward for 11 years. Initial and boundary conditions for plankton biomasses and nitrogen-based
nutrients are taken from Chenillat et al. [45]; total phytoplankton biomass and total zooplank-
ton biomass were distributed equally among initialized compartments. The ecosystem model
gets adjusted after 5 years (not shown). All the results herein are based upon the last 6 years of integration.

**Continuum size-structured plankton ecosystem model.** We use the ecosystem model of Ward et al. [24] that describes the temporal dynamics of four functional groups of phytoplankton and one zooplankton group, all further resolved into a continuum of size classes. The model and equations are specified in Ward et al. [24] and its appendix. The four phytoplankton functional types (pft) are analogs of *Prochlorococcus*, *Synechococcus*, small eukaryotes (non-diatoms), and diatoms; each pft is divided into 2 (called P1 and P2), 2 (S1 and S2), 5 (E1 to E5) and 5 (D1 to D5) size classes, respectively. Thus, 14 size classes of phytoplankton ranging from 0.5 to ~100 µm equivalent spherical diameter (ESD, (Fig 1A)) are initialized. The zooplankton are divided into eight size classes (Z1 to Z8, ranging from ~5 to 1000 µm ESD) (Fig 1A): Z1 to Z3 are nanozooplankton-like, Z4 and Z5 microzooplankton-like, and Z6 to Z8 mesozooplankton-like. Each size class is twice the volume of the previous size class. This model is based on carbon, nitrogen and iron. For their growth, phytoplankton access multiple inorganic nutrients, including nitrate (except *Prochlorococcus*), nitrite, ammonium, and iron. Model equations and parameters are taken from the Ward et al. [24] configuration that has been applied globally. We use this familiar formulation despite the many simplifications inherent to such a model structure (cf. [43]), so as to focus in detail on the structure and consequences of the grazing formulation. One distinctive feature of this model is that phytoplankton growth is a light- and temperature-dependent function of intracellular nutrient reserves, i.e., of cell quotas [48,49]. Because such a quota formulation is known to better represent plankton dynamics and nutrient cycling [50,51], the modelling community is showing growing interest in quota formulation in biogeochemical models at the global scale (e.g., with PISCES-QUOTA [18]). For clarity, below, we call the model QUOTA.

Available inorganic nutrients are used for phytoplankton growth, and phytoplankton are grazed by zooplankton. Processes such as mortality, sloppy feeding, and egestion allow the transfer of living organic material into sinking particulate detritus and dissolved organic matter. Organic detritus is returned to inorganic form through a simple parameterization of bacterial remineralization.

Following previous work that uses continuum size-structured plankton ecosystem models [21,22], the QUOTA model uses allometric relationships. Allometric relationships use a power-law function that links physiological traits *p* to cell volume *V*, thus reducing the number of free parameters in such a complex model, as follows:

$$p = aV^b \qquad (1)$$

*a* and *b* are the allometric parameters, with *b* describing the size dependency. In QUOTA, allometric relationships are used for physiological traits such as nutrient uptake, quota size, growth, mortality, sinking, and grazing rates. For most of these processes, *a* is held constant among phytoplankton groups except for maximum growth rates which vary with taxonomic group [32,52]. Within each pft, larger cells have a slower growth rate due to the allometric relationships, but among cells of similar size, diatoms will grow faster than, for example, small eukaryotes (Fig 1A). As in the original version of QUOTA, the maximum grazing rate follows allometric relationships defined by Hansen et al. [53] as a function of cell volume and we follow the Ward et al. [24] parameter choices (Fig 1A), although we find the potential maximum rates to be rather elevated in some cases. Zooplankton can feed on either phytoplankton or zooplankton, or some combination of both. Zooplankton attack rate is not modified by the density of individual prey types [39]; however, zooplankton feed preferentially on the denser prey between phytoplankton and zooplankton. Each zooplankton size class has a prey size
<figure><img src="image_4.png" /><figcaption>Fig 1. (a) Maximum phytoplankton growth rates and maximum zooplankton grazing rates as a function of cell size (volume and equivalent spherical diameter, ESD) and taxon. (b) Schematic representation of the ecosystem model simplified to show only grazing and predation fluxes for clarity (for a complete representation of the model, see Fig 2 in Ward et al. (2012) [24]). Phytoplankton functional types (Prochlorococcus (P), Synechococcus (S), small Eukaryotes (E), Diatoms (D)) are divided into different size classes (P1, P2, S1, S2, etc.). Phytoplankton and zooplankton size classes (Z1, Z2, Z3, etc.) are represented with a color pallet used in subsequent figures.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0252033.g001

preference because smaller prey are less likely to escape and are easier to ingest. A size-specific predator can thus feed on a variety of prey (zooplankton or phytoplankton from different pft, Fig 1B). This preference is set to a log-normal distribution centered to a predator to prey length scale of 10 and a standard deviation of 0.5 [11].

**Grazing refuge formulation.** The specific grazing rate *G* of a predator *jpred* on a prey *jprey* is given by:

$$G_{jpred,jprey} = G_{jpred}^{max} \cdot R \cdot \Phi_{(P \text{ or } Z)} \quad (2)$$
<table><caption>Table 1. Various mathematical formulations of the grazing rate R.</caption><thead><tr><th>Case</th><th>Grazing refuge function</th><th>R(X) mathematical formulation</th><th>parameters to set</th><th>values</th><th>units</th></tr></thead><tbody><tr><td>1</td><td>No refuge</td><td>R(X) = <sup>X</sup>&frasl;<sub>X+k<sub>pred</sub></sub></td><td>-</td><td>-</td><td>-</td></tr><tr><td>2</td><td>Ivlev (Iv)</td><td>R(X) = <sup>X</sup>&frasl;<sub>X+k<sub>pred</sub></sub> (1 - e<sup>AX</sup>)</td><td>Λ</td><td>-1</td><td>-</td></tr><tr><td>3a</td><td rowspan="3">Sigmoid (Sig)</td><td rowspan="3">R(X) = <sup>X</sup>&frasl;<sub>X+k<sub>pred</sub></sub> <sup>X<sup>2</sup></sup>&frasl;<sub>X<sup>2</sup>+K<sup>2</sup></sub></td><td rowspan="3">K<sup>2</sup></td><td>2</td><td rowspan="3">(mmol C m<sup>-3</sup>)<sup>2</sup></td></tr><tr><td>3b</td><td>4</td></tr><tr><td>3c</td><td>6</td></tr><tr><td>4a</td><td rowspan="4">Ivlev Size-Independent-Threshold (Iv-SIT)</td><td rowspan="4">R(x) = {<sup>X-T</sup>&frasl;<sub>(X-T)+k<sub>pred</sub></sub> (1 - e<sup>AX(T-T)</sup>), X &gt; T0, X &le; T</td><td rowspan="4">T</td><td>0.25</td><td rowspan="4">mmol C m<sup>-3</sup></td></tr><tr><td>4b</td><td>0.50</td></tr><tr><td>4c</td><td>0.75</td></tr><tr><td>4d</td><td>1.00</td></tr><tr><td>5a</td><td rowspan="2">Ivlev Size-Dependent-Threshold (Iv-SDT)</td><td rowspan="2">R(x) = {<sup>X-T</sup>&frasl;<sub>(X-T)+k<sub>pred</sub></sub> (1 - e<sup>AX(T-T)</sup>), X &gt; T0, X &le; T</td><td rowspan="2">a<sub>T</sub>, b<sub>T</sub></td><td>0.012, 0.27</td><td rowspan="2">mmol C m<sup>-3</sup>, -</td></tr><tr><td>5b</td><td>0.025, 0.27</td></tr></tbody></table>

For clarity, R(X) is illustrated for a single prey X given in biomass (mmol C m⁻³), although all models include the same multiple prey size classes. Bold numbers correspond to cases presented in the core paper (other cases are presented in supporting information). Λ is the rate at which saturation of the ingestion is achieved with increasing food level. T is the zooplankton feeding threshold. a_T and b_T represent allometric parameters for the size-dependent feeding threshold.

https://doi.org/10.1371/journal.pone.0252033.t001

where $G_{pred}^{max}$ is the maximum grazing rate on prey *jprey* (referenced as the ingestion rate in Ward et al. [24]), Φ is the prey biomass (which can be phytoplankton *P* or zooplankton *Z*); *R* describes the prey refuge as a function of the prey availability and reduces the grazing effort as prey biomass becomes scarce: this function typically has a sigmoidal shape. Note that although we illustrate Eq (2) in terms of a single prey item, all models include multiple prey classes and the thresholds refer to summed concentration of all suitable prey combined, weighted by the size preference function described above. In this study, we evaluate the sensitivity to the specific formulation of the prey refuge by modifying the shape function *R* as described in Table 1 and represented in Fig 2 (Note that, for clarity, Table 1 illustrates the prey refuge for a single prey. The more general formulation is given in S1 Table). We consider five different cases:

* In **case 1** (No refuge), we investigate how removing all forms of prey refuge, using only the rectangular hyperbola formulation for grazing passing through (0,0), shapes resulting ecosystem structure. In all subsequent cases, a grazing refuge is introduced by a function applied to the rectangular hyperbola in case 1.

* In **case 2** (Iv), we test the original version of the prey refuge from QUOTA that multiplies the rectangular hyperbola by an Ivlev [54] function as done by Ward et al. [24].

* In **case 3** (Sig), we replace the Ivlev reduction by a sigmoidal function as an alternative representation of a prey refuge.

* In **case 4** (Iv-SIT, for Ivlev Size-Independent Threshold), we add an explicit prey refuge to the Ivlev reduction from case 2 with a fixed, size-independent threshold prey concentration *T* [39] (see Fig 2B, horizontal lines).

* In **case 5** (Iv-SDT, for Ivlev Size-Dependent Threshold), we introduce a threshold prey concentration *T* that varies with the body size of the grazer (Fig 2B, curved lines).

It should be noted that cases 2–5 utilize variations of Ivlev or sigmoidal functions only as alternative mathematical means to reduce grazing effort at low prey concentrations and are not meant to correspond to a mechanistic description of prey encounter or feeding effort.
a

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

b

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

**Fig 2.** (a) Grazing effort as a function of prey biomass for cases 1–4 from Table 1. Solid lines correspond to cases presented in the core paper (other cases are presented in supporting information). (b) Specification of an explicit refuge food concentration T where grazing is initiated for case 4 (size-independent threshold, green lines with four different parameter values) and case 5 (body size-dependent threshold T as a function of zooplankton biovolume and ESD, pink lines). Black and white dots correspond to low and high estimates from the literature, respectively; pink lines represent T fitted to these estimates with an allometric relationship T = a_T V^b, with a_T-low = 0.012 mmol C m⁻³, a_T-high = 0.025 mmol C m⁻³, and b = 0.27 (cases 5, Iv-SDT). Green lines represent size-independent T (cases 4, Iv-SIT).

https://doi.org/10.1371/journal.pone.0252033.g002

These formulations of the prey refuge R in Eq (2) do not correspond to the classical functional response curves (e.g., [39]).

Feeding thresholds have been defined differently in the literature, depending on whether a grazer's ingestion rate or clearance rate is considered. For ingestion rates, a threshold is considered the lowest prey concentration at which positive ingestion is measured (e.g., *T*ᵢ in S1A Fig). Examination of the experimental literature shows a broad range of numerical values for such a threshold, with a mean of 1.6, minimum of 0.2 and maximum of 2.6 mmol C m⁻³ for copepods and other crustaceans, ciliates, and heterotrophic dinoflagellates [55–69].
In the context of clearance rates, a threshold is a prey concentration at which clearance rate (i.e., volume swept clear of prey predator⁻¹ time⁻¹) decreases below its maximum, as prey concentrations decline (T_c in S1B Fig). From the same literature sources as above, such thresholds show a mean of 2.4, a minimum of 0.4 and a maximum of 5.0 mmolC m⁻³. This broad range of measured values leads us to consider two alternative formulations for feeding thresholds. In the first, the grazing threshold is held constant for all size classes of zooplankton (case 4, illustrated in Fig 2B). Different values for this fixed threshold are considered, ranging from 0.25 to 1.0 mmolC m⁻³ (Fig 2A and 2B). In the second, the grazing threshold increases with increasing body size of the grazer (case 5a and 5b, Fig 2B). For case 5, we treat an example of low body size-dependence (i.e., the threshold varies from 0.04 to 2.9 mmolC m⁻³) and high body size-dependence (threshold from 0.08 to 6.0 mmolC m⁻³) across consumer biovolumes ranging from 10² to 10⁹ µm³ (Fig 2B).

In cases 3 and 4, we test the sensitivity to the prey refuge parameters. However, for simplicity, we choose to present only one of the three and four experiments performed in **cases 3 and 4**, respectively, in the core paper: case 3b and 4c (with sensitivity tests results provided in supporting information).

**Metrics for quantifying ecosystem structure and response to grazing formulation**

**Chl-a spatial distribution.** As a first index of the ecosystem model response to grazing formulation, we choose Chl-*a* annual climatological mean concentration. The model Chl-*a* is compared to satellite and in situ data. We use climatological SeaWiFS (Sea-Viewing Wide Field-of-view Sensor) data averaged over 1997–2005. We also examine the vertical and cross-shore structure of the model using CalCOFI in situ chlorophyll-*a* climatologies (https://calcofi.org/) averaged over 1984–2000, along line 70 (south of Monterey, California), from the coast to 350 km offshore. Climatologies of total phytoplankton and zooplankton carbon biomass were obtained from data reported in Kenitz et al. [70], along CalCOFI line 90 for the period Nov. 2004-Dec. 2010 for phytoplankton and Nov. 2005–Jan. 2011 for mesozooplankton. Phytoplankton were sampled at 3 depths in the euphotic zone by water bottles and analyzed by flow cytometry and epifluorescence microscopy, then each taxon converted to organic carbon from literature regressions (see [70,71] for details) and vertically integrated from the deepest sampling depth to the surface. Unlike Kenitz et al. [70], we utilize all 3 sampling depths. For microplankton taxa known to include mixotrophs, half of the biomass was assigned to phytoplankton and half to zooplankton carbon. Mesozooplankton were sampled by a vertical PRPOOS net (Planktonic Rate Processes in Oligotrophic Ocean Systems: 0.5 m diameter, 202 µm mesh) from 210–0 m, preserved in 1.8% buffered formaldehyde, Zooscanned [72,73], converted to organic carbon from body length-to-carbon regressions [74], and vertically integrated over the upper 210 m.

**Ecosystem structure and function.** We then examine the response of the planktonic ecosystem in terms of phytoplankton functional types and both phytoplankton and zooplankton size classes. In each case, we consider that a type or size class emerges if its biomass exceeds a threshold of 10⁻⁵ mmol C m⁻³ on each cell grid. We compute the Z:P carbon biomass ratio, which provides insight into the structure and function of the plankton ecosystem.

We consider the seasonal variation of biomasses and grazing fluxes in order to diagnose the emergence of community succession and understand the effect of grazing pressure on plankton ecosystem structure. The seasonal development of the different quantities (biomasses, fluxes, Z:P ratio) is computed over the full depth, excluding up to 450 km (30 grid points) at the edge of the modeled domain to remove noisy signals generated by boundary forcing.
**Spatial diversity and niches.** The model permits the coexistence of different plankton types in the same location [24]. To diagnose such coexistence, we use species richness. We define species richness as the number of phytoplankton size classes of each pft, or number of zooplankton size classes, given by the biomass concentration *Bj*, that exceeds a relative threshold of $Bj>10^{-5} \max(Bj)$ [26].

Environmental heterogeneity of physical properties creates a variety of ecological conditions that promote different habitats favoring different plankton types. In addition to these environmental factors, functional traits of species—such as resource uptake and grazing—shape different habitats or niches. From our model outputs, we defined niches as the geographical extent where the biomass of a phytoplankton size class of each pft, or zooplankton size class, dominates. Both richness and niches are computed from annual mean biomass from the 6-year averages.

# Results

## Chl-*a* spatial distribution

The annual mean chlorophyll concentration across the California Current System from SeaWiFS satellite observations (Fig 3 “Obs.”, column 1) shows the expected band of elevated Chl-*a* near the coast, sharp lateral gradients in the offshore direction, and a diminution in waters south of Pt. Conception (~34°N). The vertical section of mean Chl-*a* measurements from CalCOFI line 70 (~35.5°N, Fig 3 “Obs.”, column 2), just south of Monterey Bay, shows a near-shore surface maximum, with decreasing concentrations and the standard deep Chl-*a* maximum (DCM) in the offshore at 70–80 m depth. The five model cases with different grazing formulations differ markedly in relation to these basic geographic and vertical patterns (Fig 3). Cases 1 (*no refuge*) and 2 (*Iv*) present the largest disagreement with observations, resulting in pronounced overestimation and underestimation of chlorophyll-*a*, respectively. These biases are apparent in both the offshore and southern extent of elevated Chl-*a* concentrations seen in plan view and misrepresentation of the vertical Chl-*a* structure in section view (Fig 3.1 and 3.2, column 1 & 2). In contrast, modeled chlorophyll concentrations from grazing cases 3 (*Sig*), 4 (*Iv-SIT*) and 5 (*Iv-SDT*) show less pronounced difference in relation to observations from SeaWiFS and CalCOFI (Fig 3.3b, 3.4c, 3.5a and 3.5b, columns 1 & 2). The cross-shore gradient is somewhat better represented, with elevated concentrations at the coast extending to lower concentrations in offshore oligotrophic waters, and a DCM at about 70–80 m depth. We wish to underscore that the only variation among these different cases is altered grazing formulations.

## Ecosystem structure and seasonal cycle

Among the 22 plankton compartments initialized in the model (including both size and pft), 14 to 18 plankton compartments (i.e., 64 to 82% of initialized plankton compartments) can emerge and persist, depending on the grazing dynamics. The two largest *diatoms* (D4 and D5) and the two largest *small eukaryotes* (E4 and E5) never persist, except in case 2 where E4 faintly emerges (concentrations <10⁻⁷ mmol C m⁻³ on average; details not shown).

A seasonal cycle of biomass of varying amplitude occurs in all grazing cases, but maxima can vary 20-fold depending on the grazing formulation (compare phytoplankton concentrations in Fig 4 between cases 1 and 2; note the change in scale). In terms of phytoplankton dominance, in case 2, *Prochlorococcus* is the most abundant pft except in September, when *small eukaryotes* have a slight advantage. Note that in case 1 *diatoms* barely occur, and in cases 2–4 *diatoms* have the smallest biomass (S2 Table, S2 Fig). In contrast with case 2, total phytoplankton biomass is higher in cases 3 and 4 by a factor of 2 and 4, respectively
Obs.

1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

<figure><img src="image_21.png" /><figcaption>Chlorophyll-a concentration (mg Chl-a m⁻³)</figcaption></figure>

2

3b

4c

5a

5b
**Fig 3.** (column 1) Annual mean surface chlorophyll-a concentration (mg Chl-a m⁻³); the top panel corresponds to observations from SeaWiFS, while other panels correspond to model outputs for grazing cases 1 to 5b. (column 2) Vertical sections of annual mean chlorophyll-a concentration across CalCOFI line 70 (mg Chl-a m⁻³); the top panel corresponds to observations from CalCOFI, while other panels correspond to model outputs at the same location for cases 1 to 5b.

https://doi.org/10.1371/journal.pone.0252033.g003

(Figs 3 and 4). The relative dominances of phytoplankton pfts are similar between cases 3 and 4: all year, *Prochlorococcus* biomass is higher than *Synechococcus* biomass, followed by the biomass of *small eukaryotes* and then *diatoms*. Lastly, in both cases 5a and 5b, the dominant biomass shifts markedly to large size classes of *diatoms* and *small eukaryotes*, which have twice the biomass of *Synechococcus* and *Prochlorococcus*. Among the dominant pfts, *diatoms* can dominate during spring, while *small eukaryotes* outcompete *diatoms* during fall and winter (Fig 4).

In terms of phytoplankton size structure, for all cases and within each pft, the amplitude of biomass variation decreases inversely with size, with the two smallest size classes emerging systematically as dominants (Fig 5, S2 Table, S2 Fig). The bigger size classes of *small eukaryotes* and *diatoms* can emerge only in certain grazing pressure conditions that include explicit refugia (see below), which allow significant emergence of 5 to 8 size classes of zooplankton (Figs 4 and 5). Despite relatively small differences in emergence of phytoplankton size classes (S2 Fig), the zooplankton size structure varies greatly among cases (Fig 4). In cases 1 and 2, large zooplankton (mesozooplankton-like, Z6, Z7, or Z8) dominate throughout the year; in case 3 there is a succession of dominance between microzooplankton during spring, and larger zooplankton during summer and fall. In all subsequent cases, there is a shift to smaller zooplankton size classes. In cases 4c and 4d, the smallest zooplankton (nanozooplankton) dominate all year long (Fig 4 and S3 Fig), or alternate with mid-sized zooplankton in cases 4a and 4b (S3 Fig); in case 5, mid-sized zooplankton dominate. Independent of dominance, all size classes of zooplankton show a clear seasonal progression from spring to fall (Fig 4). For cases 2 to 5, zooplankton biomass reaches its maximum in September for large-sized zooplankton, and between June and August for mid-sized zooplankton. However, small-sized zooplankton biomass peaks at different dates depending on the grazing cases: in March in cases 2 and 3, in June in case 4 and in October in case 5.

The CCS is known to present a strong cross-shore gradient of physical and biogeochemical properties, with contrasting ecosystem structure (dominated by different sizes (and pfts) of plankton). In order to check whether the ecosystem was sensitive to the cross-shore gradient of environmental conditions, we performed similar diagnostics separately on nearshore and offshore regions for the Southern CCS (S8-S10 Figs for phytoplankton and zooplankton, respectively). For all experiments, in terms of concentrations and seasonal dynamics, we find that the ecosystem responds positively to coastal nutrient input from the coastal upwelling with earlier and higher production compared to the offshore regions. However, in both regions we find similar sensitivity to the formulation of grazing refugia as we observe for the average over the entire CCS (Fig 4).

## Trophic interactions and variability of grazing fluxes

In all cases, grazing fluxes show a clear seasonal cycle in phase with the temporal evolution of zooplankton biomass (S4 Fig), and grazing on phytoplankton dominates over predation on zooplankton by a factor of about 2. Here, we mainly discuss the spatio-temporal average of trophic interactions, as illustrated in Fig 5, which summarizes emergent compartments and grazing fluxes (see S4 Fig for detailed seasonal cycles):
<figure><img src="image_22.png" /></figure>
<caption>**Fig 4.** Seasonal variation of (column 1) the 4 plankton functional types (mmol C m⁻³) and (column 2) detailed size classes of zooplankton (mmol C m⁻³).</caption>

https://doi.org/10.1371/journal.pone.0252033.g004

* Case 1 is a peculiar case highlighting emergence of a simple trophic assemblage with only three primary compartments: two phytoplankton size classes (E1 and less importantly S2) and one zooplankton size class (Z6; see Fig 5.1). Note that Z3 emerges weakly but represents a key compartment because it involves fluxes of carbon: grazing of E1 by Z3 and predation on Z3 by Z6. Grazing of S2 by Z3 is one order of magnitude lower. It also reveals that the emergence of large zooplankton (Z6) is permitted by predation on small zooplankton (Z3).

In contrast, all other cases present a more complex response of the plankton ecosystem characterized by a large number of significant grazing and predation fluxes (more than ten):

* In case 2, four phytoplankton size classes (P1, S1, S2, E1) and two zooplankton size classes (Z6 and Z7) dominate the ecosystem. However, intermediate to high fluxes through consumers involve all the zooplankton size classes from Z1 to Z6 (dominant predation fluxes in Fig 5). Even though Z1 to Z5 do not dominate the ecosystem in terms of biomass, they play a key role in ecosystem structure. Indeed, the smallest zooplankton (Z1 to Z5) (i) grazes on the smallest size classes of *Prochlorococcus*, *Synechococcus*, and *small eukaryotes* and (ii) are predated by large zooplankton, which can then emerge.

<figure><img src="image_23.png" /><figcaption>**Fig 5.** Simplified schematic representation of grazing and predation fluxes in the ecosystem model for grazing cases 1 to 5b. For each case, biomasses (circles) and fluxes (arrows) are color-coded as follows: Dominant in red, intermediate in blue, and low in dark grey. Negligible fluxes are not shown, and negligible biomasses are represented in light grey. See S1 and S2 Tables for detailed biomasses and grazing fluxes, respectively.</figcaption></figure>

https://doi.org/10.1371/journal.pone.0252033.g005
* Case 3b shifts to a broader spectrum of co-occurring zooplankton, with dominance of zooplankton from small Z1 to large Z7, and only three phytoplankton compartments of small size (the smallest size classes of *Prochlorococcus*, *Synechococcus*, and *small eukaryotes*). The primary fluxes involve grazing of small-sized organisms (grazing of P1 by Z1, of P2 by Z2) and predation of Z1 by Z4. Even though grazing of larger phytoplankton size classes (*small eukaryotes* and *diatoms*) and predation of mid-size zooplankton are low, they are not negligible, showing that the functioning of the ecosystem involves a large number of fluxes (N = 12).

* In case 4c, emergence of larger zooplankton size classes is inhibited, limiting zooplankton diversity. Small size classes of zooplankton and phytoplankton dominate the ecosystem (Z1 to Z4, and P1, P2, S1, and E1, respectively), with a reduced number of fluxes (8).

* Case 5a contrasts with all other cases: ecosystem structure is dominated by large phytoplankton (S2, E1, E2, E3, D1, and D2) and by small and mid-sized zooplankton (Z2, Z3, and Z4). Notably, this is the first case where any size class of *diatoms* emerges as a dominant phytoplankton leading to major carbon flux. There are more dominant grazing fluxes than in other cases, and they involve grazing by zooplankton Z1 to Z5 on phytoplankton P1 to D2. Predation terms are significant and concern fluxes from Z1, Z2, to Z4, Z5, respectively. In this case, the functioning of the ecosystem is more complex than in preceding ones, with numerous grazing and predation fluxes (ten in total). Case 5b is similar to case 5a with a slight shift to larger size classes for phytoplankton (S2 is no longer dominant), significant diatom fluxes (D1 and D1), and larger zooplankton (Z2 is not dominant anymore while Z5 turns to dominant).

The relative zooplankton and phytoplankton biomass can be assessed with Z:P ratio (Fig 6) and compared with in situ observations from CalCOFI line 90 (33.485°N, 117.768°W to 30.418°N, 123.999°W) in the southern part of our study domain (Fig 6 “Observations” in top row). In the peculiar case 1, with no grazing refuge of any kind, Z:P ratio is very low (0.1) all year long and over the entire domain, suggesting weak grazing control of primary production and marked departure from observations. For cases 2 to 5, Z:P ratios show a seasonal cycle with variable relative values: in cases 2 and 3b, the Z:P ratio is high (~2) and consistently exceeds observations. In case 4c the Z:P ratio is low (= 0.8). In cases 5, the Z:P ratio is increased compared to case 4, resulting in an intermediate Z:P ratio and seasonal cycle that provides somewhat better agreement with observations. The difference between parameter values in 5a and 5b results in a higher Z:P ratio in case 5b, with an overall seasonal mean of 1.6, which compares with a value of 1.5 from observations (Fig 6). From a spatial point of view, all cases (except case 1) present a cross-shore gradient with higher Z:P at the coast, but the cases differ substantially in the spatial extent of cross-shore extension of Z:P ≥ 1, which varies from ~100 km to 1000 km offshore.

Again, we check whether the total biomasses and Z:P ratio are sensitive to the cross-shore gradient of environmental conditions (S10 Fig). We find that total phytoplankton and zooplankton biomasses averaged over the nearshore and offshore areas for the Southern CCS respond coherently to the cross-shore gradient with higher concentrations inshore in response to higher concentrations of upwelled nutrients. Apart from the extreme case 1, the nearshore and offshore regions exhibit similar sensitivities to grazing dynamics.

## Spatial niches and richness

Both total species richness and the spatial niches occupied by different plankton categories are highly sensitive to the zooplankton grazing formulation used. In the peculiar case 1, in terms
Observations
CalCOFI line 90

(m m-3) ssaamqoq mmm

1

2

3b

4c

5a

5b

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

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

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

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

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

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

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

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

Z:P ratio

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

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

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

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

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

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

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

Z:P RATIO
**Fig 6. Z:P ratio.** (column 1) Seasonal variation of total phytoplankton biomass (blue), total zooplankton biomass (red) and the associated Z:P ratio (black), in mmolC m⁻³. The top panel corresponds to observations from CalCOFI line 90, while other panels correspond to model outputs for grazing cases 1 to 5b. The number in the upper-left corner indicates the seasonal mean of the Z:P ratio. (column 2) Spatial distribution of the annual-mean Z:P ratio. The solid black line indicates Z:P = 1.

https://doi.org/10.1371/journal.pone.0252033.g006

of equilibrium richness, up to 9 phytoplankton compartments and 8 zooplankton compartments emerge (Fig 7, top row of columns 1 & 2), while niches are overwhelmingly dominated by only E1 and Z6 over the entire CCS (Fig 7, top row of columns 3 & 4).

In all other cases, richness decreases in the cross-shore direction, with no clear latitudinal trend (Fig 7, column 1 & 2). Species richness varies among cases, with highest richness for case 2 for both phytoplankton and zooplankton, lowest phytoplankton richness in case 4, and lowest zooplankton richness in case 5. Except in case 1, niches show a cross-shore gradient, with a coastal dominance of larger size classes for both phytoplankton and zooplankton (Fig 7, column 3 & 4). In case 2, small species niches dominate the domain with E1 in coastal waters (exclusively in the central CCS), and P1 elsewhere; despite few phytoplankton niches, there is a high diversity of zooplankton niches with large size classes (Z6 and Z7) at the coast and small ones offshore (Z1 and Z2). Case 3 presents a single dominant phytoplankton niche (P1). Despite low diversity of phytoplankton niches in both cases 2 and 3b, these cases present a high diversity of zooplankton niches (Z5 and Z7 at the coast, Z1 and Z4 offshore). In case 4, niche diversity of phytoplankton slightly increases, while niche diversity of zooplankton is the lowest among all cases. In contrast, the simulation presenting the most diverse type of phytoplankton niches is case 5a, with—from coast to offshore water—diatoms (D2 close to the coast and D1), small eukaryotes (E1), *Prochlorococcus* (P1) and *Synechococcus* (S1). Notably, not until case 5, which incorporates a body-size dependent grazing threshold, do diatoms dominate in any part of the CCS. Despite such diversity of phytoplankton niches, zooplankton niche diversity is moderate, with niches of small zooplankton (nano and micro) dominating at the coast (niches of Z4, Z2, and Z1 from the coast to offshore).

Discussion

This analysis illustrates that a suite of properties relevant to ocean biogeochemistry, pelagic food web structure, and carbon flow pathways are highly sensitive to the formulation of grazing in a contemporary planktonic ecosystem model. The model structure incorporates multiple forms of dissolved nutrients, intracellular nutrient quotas, four phytoplankton functional types, and a spectrum of size classes of both phyto- and zooplankton, with allometrically scaled parameters. The food web model is coupled to a regional ocean model (ROMS-CROCO) of the California Current coastal upwelling ecosystem, where spatial gradients are high in both the cross-shore and vertical dimensions. The model results show that changes in the functional form of grazing that are all within the scope of previous measurements can have substantial consequences for spatial patterns of plankton concentrations, species richness, food web pathways, and underlying seasonal cycles.

Our results are not unexpected, in light of the previous numerical experiments [41,42, and others]. Our findings also reinforce the results of Anderson et al. [27], who in a global ocean circulation model found high sensitivity of phytoplankton functional group distributions and biomass to the choice of grazer functional responses. However, that study did not specifically explore the importance of prey refugia or vertical distributions and considered only two zooplankton classes.

Here we discuss the results of our sensitivity analyses with increasing complexity of the grazing formulation, from ignoring refugia to implementation of a size-dependent grazing
RICHNESS

1

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

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

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

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

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

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

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

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

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

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

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

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

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

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

NICHES

Phytoplankton

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

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

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

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

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

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

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

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

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

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

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

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

2

3b

4c

5a

5b

<caption>Fig 7. (columns 1 and 2) Richness and (column 3 and 4) niches of phytoplankton (columns 1 and 3) and zooplankton (columns 2 and 4) for grazing in cases 1 to 5b.</caption>

https://doi.org/10.1371/journal.pone.0252033.g007

refuge. Our principal objective was not to tune parameters to most closely resemble observations but to understand the sensitivity to the form of grazing considered and provide guidance on preferred formulations of grazing in the future.

**No refuge = no diversity**

*Case 1* eliminates all forms of grazing refuge, replacing it with a classic rectangular hyperbola that goes through the origin (0, 0). In *case 1*, model output disagreed with observations in the
CCS with a substantial biomass overestimation, almost 5 times over chlorophyll-a observations. Not surprisingly, the formulation without grazing refuge fails to reproduce coherent ecosystem structure in our study region: the dominance of a single phyto- and zooplankton species and the prevalence of single species niche over the entire CCS differ completely from observations in such an upwelling system, with a cross-shore shift from coastal rich water to offshore oligotrophic waters typical of coastal upwelling systems [70,75]. Observing such a homogeneous pattern of niches in such a heterogeneous environment [45] highlights the importance of more realistic grazing formulations. Persistently low values of the Z:P ratio (spatially and seasonally) indicate abnormally high phytoplankton biomass over zooplankton biomass. Independently of the degree of realism of ecosystem dynamics of the CCS, withdrawing a grazing refuge leads to the emergence of a short food web that involves one of the smallest phytoplankton (E1) feeding one type of nanozooplankton (Z3) which is predated by a single type of mesozooplankton (Z6). The most competitive phytoplankter is not the smallest—with the highest nutrient affinity - or the biggest—with the highest growth rate—but an intermediate. In sum, in this peculiar case 1, by removing the grazing refuge, we introduce non-negligible grazing pressure over low-density prey that limits the emergence of an equal diversity among plankton size classes independently of environmental variability. In the HNLC (High Nutrient-Low Chlorophyll) conditions in the Equatorial Pacific, Leising et al [41] similarly could not reproduce ambient conditions without inclusion of a grazing threshold in a simpler Fe-N-P-Z model.

## Classic refuge formulations and the emergence of diversity

Empirical formulations of grazing that are typically used in modelling studies include a simple and fixed grazing refuge for prey (see [42] for references, or as in the original version of the QUOTA model from [24]). We test the sensitivity of such empirical formulation in *case 2* and *case 3*. Implementing such a refuge for low prey density with classic sigmoidal grazing formulation increases phytoplankton diversity of smaller species, i.e., *Prochlorococcus* and *Synechococcus*, as well as zooplankton diversity. This diversity mirrors complex trophic interactions with multiplicity of grazing and predation fluxes from small to large size classes. Interestingly, diversity of zooplankton is dominated by two large zooplankton in case 2, while in case 3, this diversity is homogeneous with clear successions of diverse zooplankton assemblages (small and large size classes). Such an emergence of large size classes of zooplankton is possible through high predation pressure. Smaller-sized zooplankton do not dominate in case 2 because of high fluxes through this compartment (from grazers to their predators). As discussed in Ward et al. [24], a model that explicitly resolves the grazer community allows the “kill-the-winner” mechanism to maintain phytoplankton diversity, according to Armstrong [76]. However, our results reveal that such diversity depends not only on the diversity of the grazers but also crucially on how they graze by introducing a classic refuge formulation.

Such classical grazing formulations allow the emergence of diversity with limitations. With an Ivlev formulation (*case 2*), biomasses are markedly underestimated compared to observations, and larger species such as diatoms do not emerge in nutrient-rich coastal water [77]. On the other hand, with a stronger suppression of grazing at low concentrations achieved by a sigmoidal formulation (*case 3*), the misfit to phytoplankton biomass and its spatial distribution is less extreme, though still highly biased. In this case, trophic interactions are modified, leaving more capacity for small size classes compared to mid classes to grow (*case 2*). The main differences between these two cases occur within zooplankton dynamics, and not within phytoplankton dynamics, with a partial reduction in predation pressure for the sigmoid formulation.
Within *case 3*, increasing the refuge (with increased values of $K^2$) allows a reduction of grazing pressure for the same prey density. This promotes an increase of both phytoplankton and zooplankton biomass, especially emergence of small phytoplankton and nanozooplankton earlier in the year and a decrease of the Z:P ratio.

Simple refuge formulations as in cases 2 and 3 are typically used because of the lack of data [42], while it is now clear that such a grazing formulation in a continuum size-structured ecosystem model leads to substantial inaccuracies.

## Toward an explicit threshold

Lastly, we implement an explicit threshold prey concentration below which grazing is zero, as either a fixed value for all size classes (*case 4*) or a refuge that increases as the consumer size increases (*case 5*). In both *cases 4* and 5, trophic interactions are still quite complex with various emergences of pfts and size classes. We recognize that Leising et al. [41], as expanded on more fully in Strom et al. [78], illustrate how the feeding behavior of individual grazer types, each with different feeding thresholds, can appear in the aggregate to exhibit no threshold behavior.

A body size-independent grazing refuge threshold (*case 4*) increases phytoplankton and zooplankton biomass to somewhat more realistic levels and shifts plankton composition towards smaller size classes compared to other cases; it also allows an increase in phytoplankton niches and a decrease in zooplankton niches. Increasing the threshold refuge *T* (*cases 4*) enhances phytoplankton biomass (*S5* and *S6 Figs*) through less top-down control (i.e., less grazing pressure through less zooplankton biomass, (*S6 Fig*)) but decreases both niche diversity and richness for phyto- and zooplankton (*S7 Fig*).

Implementing a size-dependent refuge threshold (*case 5*) allows significant changes in plankton dynamics: (i) there is a somewhat better agreement of phytoplankton biomass with observations and (ii) plankton ecosystem structure is more realistic for the CCS with the emergence of large phytoplankton and, notably, coastal dominance of diatoms (for *case 5a*) [77]. In this case, the larger the plankton size class is, the higher the refuge. Zooplankton diversity is relatively low and is mainly made up of smaller zooplankton (nano- and microzooplankton), involving more grazing fluxes than predation fluxes compared to other cases. There are more trophic interactions (or fluxes) at play but with lower intensity. Altered parameter values within case 5 mainly affect phytoplankton, permitting larger organisms to emerge. The steeper body size-dependent thresholds enable both small and large organisms to escape grazing and finally allow diatoms to dominate in the ecosystem, promoting diversity (i.e., in terms of richness and diversity of niches). This leads to higher carbon transfer into the food web as soon as the bloom of diatoms appears (higher Z:P ratio), and permits more realistic Z:P ratios in relation to in situ observations. Some evidence for a body size-dependent feeding threshold within a species comes from the results of Olivares et al. [79], whose study suggests that adult female *Paracartia grani* have appreciably higher threshold concentrations than younger developmental stages.

In another set of experiments, we tested the sensitivity of the model to the grazing kernel shape (results not shown). Increasing the range of prey up to reasonable values (i.e., increasing the normal distribution width from 0.5 to 0.915) or changing the normal distribution to an asymmetric distribution allowing grazing on smaller prey does not appreciably change our results.

We intentionally do not consider a number of processes that could be explored further because of our focus on changing only one aspect of an existing model structure. We do not consider the effects of acclimation to prey (e.g., [12,42]), variable satiation responses,
departures of predator:prey size ratios from 10:1 [80], or differences in half-saturation constants for different prey types or different predators. We also recognize that, in addition to consumer body size, differences among taxa and differences in ocean habitat can affect prey thresholds. For example, Paffenhöfer and Stearns [81] suggested that copepods occurring in offshore and oceanic environments show feeding thresholds at lower food concentrations than more coastal copepods. Lessard and Murrell [82] indicated that the grazing threshold for microzooplankton in the oligotrophic waters of the Sargasso Sea corresponds to the annual minimum chlorophyll concentration found in the mixed layer. There are also suggestions that the minimum food concentration for positive net growth increases with temperature (e.g., [83]), although feeding thresholds may not be temperature-dependent [66]. We also do not consider the relative motility of different predator and prey types, turbulence-mediated modifications to encounter, prey switching responses, algal toxicity, or many other important details of planktonic predator-prey interactions [43,84]. Furthermore, the present model does not address the importance of the mathematical formulation of zooplankton mortality, a closure term that has been widely recognized to affect the dynamics of planktonic ecosystem models [42,85,86].

Our work does not seek to fully model realistic ecosystem structure and diversity over the CCS. Goebel et al. [87,88] utilized a "self-organizing" model [4] to describe diversity patterns off the California coast, although there were limited validation data available for the region. While the model we employed does include some diversity, this is only at a low level (only 14 phytoplankton and 8 zooplankton compartments) due to several constraints: the multiplicity of experiments conducted, the horizontal resolution of the model, and the available computing time. Nevertheless, such a simple size-structured ecosystem model allowed us to diagnose key processes.

## Conclusions

In conclusion, small changes in grazing refuge formulation drive large modifications in plankton diversity, structure, and ecosystem functioning. These results reiterate the importance of careful formulation of the grazing refuge to accurately represent plankton dynamics, as well as for a spectrum of ecosystem processes. This work seeks to inform future formulations of ecosystem models. While finer spatial resolution and better representation of mesoscale eddies and submesoscale filaments [46,89] are laudable goals for future models, here we emphasize the importance of biological structure and, in particular grazing processes. We show that fundamental differences in model results emerge from a simple modification of grazing refugia across a range of plausible functional forms.

Our results suggest that an explicit grazing refuge, with refuge prey concentration dependent on the body size of grazers, is likely to provide more coherent plankton ecosystem dynamics compared to classic formulations or size-independent threshold refugia. Therefore, our recommendation for future models of plankton ecosystem structure in a continuum size-structured model is to implement a size-dependent threshold refuge. As for parameter choice, we suggest use of our less steep size-dependent case (case 5a, i.e., $a_T = 0.012 \text{ mmol C m}^{-3}$), which appears closer to the modest available literature data. We also recommend initiation of a new suite of experimental grazing studies intended to test the hypothesis of body size-dependent grazing refugia across a broad spectrum of planktonic grazers, and to rigorously establish appropriate parameter values for future models.

One may wonder whether such sensitivity to grazing pressure would emerge in another region. In light of the environmental heterogeneity observed in the California Current Ecosystem, which extends from nutrient-rich coastal waters to offshore oligotrophic waters and
encompasses most of the range of optical properties found in the world ocean [90], we suggest that these findings are likely quite general. Moreover, given the importance of Eastern Boundary Current Upwelling Ecosystems for marine biogeochemistry and trophic interactions, it is particularly important that they be more accurately represented in models that seek to represent global scale ocean processes.

**Supporting information**

S1 Fig. Schematic representation of (a) ingestion rate and (b) clearance rate for a classic sigmoidal functional response (solid line) vs. a response considering a feeding threshold $T_i$ and $T_c$, (for ingestion rate and for clearance rate, respectively) (dashed line).
(DOCX)

S2 Fig. Seasonal variation of detailed pfts accounting for emergent phytoplankton size classes (in mmol C m⁻³).
(DOCX)

S3 Fig. Seasonal variation of (column 1) the 4 plankton functional types (mmol C m⁻³) and (column 2) detailed size classes of zooplankton (mmol C m⁻³).
(DOCX)

S4 Fig. Seasonal variation of grazing fluxes (averaged on the last year of simulation, integrated over depths and averaged over longitude and latitude, excluding boundary conditions). Grazing fluxes correspond to equation 2. Color code corresponds to each zooplankton size class (as in Fig 3) while dotted-dashed, dashed and solid lines correspond to preys of increasing maximum growth and grazing rates, respectively (y-axis in mmol C m⁻³ d⁻¹).
(DOCX)

S5 Fig. (column 1) Annual mean surface chlorophyll-a concentration (mg Chl-a m⁻³); (column 2) Vertical sections of annual mean chlorophyll-a concentration across CalCOFI line 70 (mg Chl-a m⁻³). Results correspond to model outputs for grazing cases 3 and 4 not shown in the core paper.
(DOCX)

S6 Fig. Z:P ratio. (column 1) Seasonal variation of total phytoplankton biomass (blue), total zooplankton biomass (red) and the associated Z:P ratio (black). The number in the upper-left corner indicates the seasonal mean of the Z:P ratio. (column 2) Spatial distribution of the annual-mean Z:P ratio. The solid black line indicates Z:P = 1. Results correspond to model outputs for grazing cases 3 and 4 not shown in the core paper.
(DOCX)

S7 Fig. (columns 1 and 2) Richness and (column 3 and 4) niches of phytoplankton (columns 1 and 3) and zooplankton (columns 2 and 4) for grazing in cases 3 and 4.
(DOCX)

S8 Fig. Comparison of phytoplankton concentrations between coastal (0-100km) and offshore region (300-400km), in the Southern CCS (+/- 0.5° around line 90 of CalCOFI): Seasonal variation of detailed pfts accounting for emergent phytoplankton size classes (in mmol C m⁻³).
(DOCX)

S9 Fig. Comparison of zooplankton concentrations between coastal (0-100km) and offshore region (300-400km), in the Southern CCS (+/- 0.5° around line 90 of CalCOFI):
Seasonal variation of detailed size classes of zooplankton (mmol C m⁻³).

(DOCX)

S10 Fig. Comparison of Z:P ratio between coastal (0-100km) and offshore region (300-400km), in the Southern CCS (+/- 0.5° around line 90 of CalCOFI): Seasonal variation of total phytoplankton biomass (blue), total zooplankton biomass (red) and the associated Z: P ratio (black) from observations (CalCOFI data) and the different numerical experiments. The number in the upper-left corner indicat.

(DOCX)

S1 Table. Various mathematical formulations of the grazing rate R, in their general formulation, i.e., for several prey (in biomass, mmolC m⁻³). All models include the same multiple prey size classes. Bold numbers in the first column indicate cases presented and discussed in the core paper, other cases are presented in supporting information. Λ is the rate at which saturation of the ingestion is achieved with increasing food level. T is the zooplankton feeding threshold. aT and bT represent allometric parameters for the size-dependant feeding threshold.

(DOCX)

S2 Table. Biomasses (mmolC m⁻³). Biomasses are averaged over the 6 last years of the numerical experiments and spatially computed over the full depth, excluding up to 450 km (30 grid points) at the edge of the modelled domain to remove noisy signals generated by boundary forcing. If biomass is below a threshold of 10⁻⁵ mmol C m⁻³ on each cell for a given plankton type, we assume that the type does not emerge (no number is given in the table below). Bold cases represent cases showed in the core paper.

(DOCX)

S3 Table. Grazing fluxes (mmolC m⁻³ d⁻¹). Grazing fluxes are averaged over the 6 last years of the numerical experiments and spatially computed over the full depth, excluding up to 450 km (30 grid points) at the edge of the modelled domain to remove noisy signals generated by boundary forcing. If a plankton type (prey or predator) does not emerge (see S1 Table), there is not available flux. Bold cases represent cases showed in the core paper.

(DOCX)

**Acknowledgments**

We thank B. A. Ward, who provided support during the model implementation phase. We also thank T. Gorgues, M. Sourisseau and P. Pondaven for their helpful discussions. We also thank M.R. Landry for use of pico- and microplankton data. Numerical simulations were performed on Ada of the IDRIS center. Model sources (ROMS-QUOTA) used in this manuscript are available on figshare with the identifier https://doi.org/10.6084/m9.figshare.12696422.v1.

**Author Contributions**

**Conceptualization:** Fanny Chenillat, Pascal Rivière, Mark D. Ohman.

**Formal analysis:** Fanny Chenillat.

**Funding acquisition:** Fanny Chenillat, Mark D. Ohman.

**Investigation:** Fanny Chenillat, Mark D. Ohman.

**Methodology:** Fanny Chenillat, Pascal Rivière, Mark D. Ohman.

**Project administration:** Fanny Chenillat.
**Resources:** Pascal Rivière.

**Software:** Fanny Chenillat.

**Supervision:** Pascal Rivière.

**Validation:** Fanny Chenillat, Mark D. Ohman.

**Visualization:** Fanny Chenillat.

**Writing – original draft:** Fanny Chenillat.

**Writing – review & editing:** Pascal Rivière, Mark D. Ohman.

## References

1. Boyd PW, Claustre H, Lévy M, Siegel DA, Weber T. Multi-faceted particle pumps drive carbon sequestration in the ocean. Nature. 2019; 568: 327–335. https://doi.org/10.1038/s41586-019-1098-2 PMID: 30996317

2. Runge JA. Should we expect a relationship between primary production and fisheries? The role of copepod dynamics as a filter of trophic variability. Hydrobiologia. 1988; 167/168: 61–71.

3. Lévy M., Franks PJS, Smith KS. The role of submesoscale currents in structuring marine ecosystems. Nature Communications. 2018; 9. https://doi.org/10.1038/s41467-018-07059-3 PMID: 30420651

4. Follows M, Dutkiewicz S, Grant S, Chisholm S. Emergent biogeography of microbial communities in a model ocean. Science. 2007; 315, 1843–1846. https://doi.org/10.1126/science.1138544 PMID: 17395828

5. Dutkiewicz S, Cermeno P, Jahn O, Follows MJ, Hickman AE, Taniguichi DAA, et al. Dimensions of marine phytoplankton diversity. Biogeosciences. 2020; 17: 609–634. https://doi.org/10.5194/bg-17-609-2020

6. Frouin R, Ramon D, Boss E, Jolivet D, Compiègne M, Tan J et al. Satellite Radiation Products for Ocean Biology and Biogeochemistry: Needs, State-of-the-Art, Gaps, Development Priorities, and Opportunities. Frontiers in Marine Science. 2018; 5. https://doi.org/10.3389/fmars.2018.00003

7. Houghton IA, Koseff JF, Monismith SG, Dabiri JO. Vertically migrating swimmers generate aggregation-scale eddies in a stratified column. Nature. 2018; 556: 497–500. https://doi.org/10.1038/s41586-018-0044-z PMID: 29670284

8. Humes AG. How many copepods?, p. 1–7. In Ferrari F. D. and Bradley B. P. [eds.]. Kluwer. 1994.

9. Finlay BJ, Clarke KJ. Ubiquitous dispersal of microbial species. Nature. 1999; 400:828.

10. Stoecker DK, Hansen PJ, Caron DA, Mitra A, Annual R. Mixotrophy in the Marine Plankton. 2017; p. 311–335. Annual Review of Marine Sciences, Vol 9. https://doi.org/10.1146/annurev-marine-010816-060617 PMID: 27483121

11. Kiørboe T. A Mechanistic Approach to Plankton Ecology. Princeton University Press, Princeton, NJ. 2008. https://doi.org/10.1007/s00442-007-0893-x PMID: 17990002

12. Franks PJS, Wroblewski JS, Flierl GR. Behaviour of a simple plankton model with food-level acclimation by herbivores. Mar. Biol. 1986; 91: 121–129.

13. Edwards CA, Powell TA, Batchelder HP. The stability of an NPZ model subject to realistic levels of vertical mixing. J. Mar. Res. 2000; 58, 37–60.

14. Franks PJS. NPZ Models of Plankton Dynamics: Their Construction, Coupling to Physics, and Application. Journal of Oceanography. 2002; 58, 379–38. https://doi.org/10.1023/A:1015874028196.

15. Moore JK, Doney S, Kleypas J, Glover D, Fung F. An intermediate complexity marine ecosystem model for the global domain, Deep-Sea Res. II. 2002; 49, 403–462.

16. Litchman E, Klausmeier CA, Miller JR, Schofield OM, Falkowski P G. Multi-nutrient, multi-group model of present and future oceanic phytoplankton communities, Biogeosciences. 2006; 3, 585–606. https://doi.org/10.5194/bg-3-585-2006.

17. Le Quéré CL, Harrison SP, Colin Prentice I, Buitenhuis ET, Aumont O, Bopp L, et al. Ecosystem dynamics based on plankton functional types for global ocean biogeochemistry models. Global Change Biology. 2005; 11: 2016–2040. https://doi.org/10.1111/j.1365-2486.2005.1004.x

18. Aumont O, Ethe C, Tagliabue A, Bopp L, Gehlen M. PISCES-v2: an ocean biogeochemical model for carbon and ecosystem studies. Geoscientific Model Development. 2015; 8 (8), 2465–2513. ISSN 1991-959X.
19. Kwiatkowski L, Aumont O, Bopp L, Ciais P. The impact of variable phytoplankton stoichiometry on pro- jections of primary production, food quality, and carbon uptake in the global ocean. *Global Biogeochem- ical Cycles*. 2018; 32, 516–528. https://doi.org/10.1002/2017GB005799.

20. Martini S, Larras E, Boyé A, Faure E, Aberle N, Archambault P, et al. Functional trait-based approaches as a common framework for aquatic ecologists. Limnology and Oceanography. 2020. https://doi.org/10. 1002/lno.11655

21. Moloney CL, Field JG. The size-based dynamics of plankton food webs. I. A simulation model of carbon and nitrogen flows. Journal of Plankton Research. 1991; 13: 1003–1038. https://doi.org/10.1093/ plankt/13.5.1003

22. Baird ME, Suthers IM. A size-resolved pelagic ecosystem model. Ecological Modelling. 2007; 203: 185–203. https://doi.org/10.1016/j.ecolmodel.2006.11.025

23. Banas NS. Adding complex trophic interactions to a size-spectral plankton model: Emergent diversity patterns and limits on predictability. Ecological Modelling. 2011; 222: 2663–2675. http://doi.org/10. 1016/j.ecolmodel.2011.05.018.

24. Ward BA, Dutkiewicz S, Jahn O, Follows MJ. A size-structured food-web model for the global ocean. Limnology and Oceanography. 2012; 57: 1877–1891. https://doi.org/10.4319/lo.2012.57.6.1877

25. Dutkiewicz S, Boyd PW, Riebesell U, Exploring biogeochemical and ecological redundancy in phyto- plankton communities in the global ocean. Glob Change Biol. 2021; https://doi.org/10.1111/gcb.15493 PMID: 33342048

26. Lévy M, Jahn O, Dutkiewicz S, Follows MJ, d'Ovidio F. The dynamical landscape of marine phytoplank- ton diversity. J. R. Soc. Interface. 2015; 12: 0481. https://doi.org/10.1098/rsif.2015.0481 PMID: 26400196

27. Anderson TR, Gentleman WC, and Sinha B. Influence of grazing formulations on the emergent proper- ties of a complex ecosystem model in a global ocean general circulation model. Progress in Oceanogra- phy. 2010; 87: 201–213. https://doi.org/10.1016/j.pocean.2010.06.003

28. Martz T, Send U, Ohman MD, Takeshita Y, Bresnahan P., Kim HJ, et al. Dynamic variability of biogeo- chemical ratios in the Southern California Current System. Geophysical Research Letters. 2014; 41: 2496–2501, https://doi.org/10.1002/2014gl059332

29. Chabert P, d'Ovidio F, Echevin V, Stukel MR, Ohman MD. Cross-shore flow and implications for Carbon Export in the California Current Ecosystem: a Lagrangian analysis. JGR Oceans. 2021. https://doi.org/ 10.1029/2020JC016308 PMID: 33842183

30. Pauly D, Christensen V. Primary production required to sustain global fisheries. Nature. 1995, 374: 255–257.

31. Capet X, McWilliams JC, Molemaker MJ, Shchepetkin AF. Mesoscale to Submesoscale Transition in the California Current System. Part I: Flow Structure, Eddy Flux, and Observational Tests. *J. Phys. Oceanogr.* 2008; 38, 29–43. https://doi.org/10.1175/2007JPO3671.1.

32. Litchman E, Klausmeier CA, Schofield OM, Falkowski PG. The role of functional traits and trade-offs in structuring phytoplankton communities: scaling from cellular to ecosystem level. Ecology Letters. 2007; 10: 1170–1181. https://doi.org/10.1111/j.1461-0248.2007.01117.x PMID: 17927770

33. Gao KS, Helbling EW, Hader DP, Hutchins DA. Responses of marine primary producers to interactions between ocean acidification, solar radiation, and warming. Marine Ecology Progress Series. 2012; 470: 167–189. https://doi.org/10.3354/meps10043

34. Calbet A, Landry MR. Phytoplankton growth, microzooplankton grazing, and carbon cycling in marine systems. Limnology and Oceanography. 2004; 49: 51–57.

35. Steele JH. The Structure of Marine Ecosystems. Harvard University Press, Cambridge, Massachu- setts. 1974.

36. Evans GT. Functional response and stability. Am. Nat. 1977; 111: 799–802.

37. Frost BW, Franzen N. Grazing and iron limitation in the control of phytoplankton stock and nutrient con- centration: a chemostat analogue of the Pacific equatorial upwelling zone. Marine Ecology Progress Series. 1992; 83: 291–303.

38. Myerscough MR, Darwen MJ, Hogarth WL. Stability, persistence and structural stability in a classical predator-prey model. Ecol. Model. 1996; 89: 31–42.

39. Gentleman W, Leising A, Frost B, Strom S, Murray J. Functional responses for zooplankton feeding on multiple resources: A review of assumptions and biological dynamics. Deep-Sea Research Part II: Topi- cal Studies in Oceanography. 2003; 50: 2847–2875. https://doi.org/10.1016/j.dsr2.2003.07.001

40. Fussman GF, Blasius B. Community response to enrichment is highly sensitive to model structure. Biol. Lett. 2005; 1: 9–12. https://doi.org/10.1098/rsbl.2004.0246 PMID: 17148115
41. Leising AW, Gentleman WC, Frost BW. The threshold feeding response of microzooplankton within Pacific high-nitrate low-chlorophyll ecosystem models under steady and variable iron input. Deep-Sea Research Part II-Topical Studies in Oceanography. 2003. 50: 2877-2894. https://doi.org/10.1016/j.dsr2.2003.07.002

42. Gentleman WC, Neuheimer AB. Functional responses and ecosystem dynamics: how clearance rates explain the influence of satiation, food-limitation and acclimation. Journal Of Plankton Research. 2008; 30: 1215-1231.

43. Flynn KJ, Mitra A. Why plankton modelers should reconsider using rectangular hyperbolic (Michaelis-Menten, Monod) descriptions of predator-prey interactions. Frontiers in Marine Science. 2016; 3. https://doi.org/10.3389/fmars.2016.00165

44. Shchepetkin AF, McWilliams JC. The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography- following-coordinate oceanic model. Ocean Model. 2005; 9:347-404. https://doi.org/10.1016/j.ocemod.2004.08.002

45. Chenillat F, Rivière P, Capet X, Franks PJS, Blanke B. California coastal upwelling onset variability: Cross-shore and bottom-up propagation in the planktonic ecosystem. PLoS One. 2013; 8: e62281. https://doi.org/10.1371/journal.pone.0062281 PMID: 23690935

46. Chenillat F, Franks PJS, Combes V. Biogeochemical properties of eddies in the California Current System. Geophys Res Lett. 2016; 43: 5812-5820. https://doi.org/10.1002/2016GL068945.

47. Chenillat F, Franks PJS, Capet X, Rivière P, Grima N, Blanke B, et al. Eddy properties in the Southern California Current System. Ocean Dynamics. 2018; 68: 761-777.

48. Droop MR. Vitamin B12 and marine ecology. IV. The kinetics of uptake, growth and inhibition in Monochrysis lutheri. J. Mar. Biol. Assoc. UK. 1968; 48: 689-733.

49. Geider RJ, Macintyre HL, Kana TM. A dynamic regulatory model of phytoacclimation to light, nutrients and temperature. Limnol. Oceanogr. 1998; 43: 679-694. https://doi.org/10.4319/lo.1998.43.4.0679

50. Flynn KJ. The importance of the form of the quota curve and control of non-limiting nutrient transport in phytoplankton models, Journal of Plankton Research. 2008; 30(4): 423-438. https://doi.org/10.1093/plankt/fbn007.

51. Ayata SD, Lévy L, Aumont O, Sciandra A, Sainte-Marie J, et al. Phytoplankton growth formulation in marine ecosystem models: should we take into account photo-acclimation and variable stoichiometry in oligotrophic areas?. Journal of Marine Systems, Elsevier, 2013, Journal of Marine Systems. 2013; 125: 29-40. https://doi.org/10.1016/j.jmarsys.2012.12.010.

52. Raven JA, Geider RJ. Temperature and algal growth. New Phytologist. 1988; 110:441-461. https://doi.org/10.1111/j.1469-8137.1988.tb00282.x

53. Hansen PJ, Bjornsen PK, Hansen BW. Zooplankton grazing and growth: Scaling within the 2-2,000-m microzooplankton body size range. Limnology and Oceanography. 1997; 42: 687-704. https://doi.org/10.4319/lo.1997.42.4.0687

54. Ivlev VS. Experimental ecology of the feeding of fishes. Pischepromizdat. 1955.

55. Frost BW. Effects of size and concentration of food particles on the feeding behavior of the marine planktonic copepod *Calanus pacificus*. Limnology and Oceanography. 1972; 17: 805-815.

56. Frost BW. A threshold feeding behavior in *Calanus pacificus*. Limnology and Oceanography. 1975; 20: 263-266.

57. Fernandez F. Particle selection in the nauplius of *Calanus pacificus*. Journal of Plankton Research. 1979; 1: 313-328.

58. Capriulo GM. Feeding of field collected tintinnid micro-zooplankton on natural food. Marine Biology. 1982; 71: 73-86.

59. Frost BW, Landry MR, Hassett RP. Feeding behavior of large calanoid copepods *Neocalanus cristatus* and *N. plumchrus* from the subarctic Pacific Ocean. Deep-Sea Res. 1983; 30: 1-13.

60. Ohman MD. Omnivory by *Euphausia pacifica*: the role of copepod prey. Marine Ecology Progress Series. 1984; 19: 125-131.

61. Caron DA, Goldman JC, Andersen OK, Dennett MR. Nutrient cycling in a microflagellate food chain: II. Population dynamics and carbon cycling. Marine Ecology Progress Series. 1985; 24: 243-254.

62. Kiørboe T, Mohlenberg F, Hamburger K. Bioenergetics of the planktonic copepod *Acartia tonsa*: relation between feeding, egg production and respiration, and composition of specific dynamic action. Marine Ecology Progress Series. 1985; 26: 85-97. https://doi.org/10.3354/meps026085

63. Strom SL. Growth and grazing rates of the herbivorous dinoflagellate *Gymnodinium* sp. from the open subarctic Pacific Ocean. Marine Ecology Progress Series. 1991; 78: 103-113. https://doi.org/10.3354/meps078103
64. Durbin AG, Durbin EG. Seasonal changes in size frequency distribution and estimated age in the marine copepod *Acartia tonsa* during a winter-spring diatom bloom in Narragansett Bay. Limnology and Oceanography. 1992; 37: 379-392.

65. Wlodarczyk E, Durbin AG, Durbin EG. Effect of temperature on lower feeding thresholds, gut evacuation rate, and diel feeding behavior in the copepod *Acartia hudsonica*. Marine Ecology Progress Series. 1992; 85: 93-106. https://doi.org/10.3354/meps085093

66. Pahlow M, Prowe AEF. Model of optimal current feeding in zooplankton. Marine Ecology Progress Series. 2010; 403: 129-144. https://doi.org/10.3354/meps08466

67. Saage A, Vadstein O, Sommer U. Feeding behaviour of adult Centropages hamatus (Copepoda, Calanoida): Functional response and selective feeding experiments. Journal of Sea Research. 2009; 62: 16-21. https://doi.org/10.1016/j.seares.2009.01.002

68. Takahashi K, Tiselius P. Ontogenetic change of foraging behaviour during copepodite development of Acartia clausi. Marine Ecology Progress Series. 2005; 303: 213-223. https://doi.org/10.3354/meps303213

69. Almeda R, Augustín CB, Alcaraz M, Calbet A, Saiz E. Feeding rates and gross growth efficiencies of larval developmental stages of Oithona davisae (Copepoda, Cyclopoida). Journal of Experimental Marine Biology and Ecology. 2010; 387: 24-35. https://doi.org/10.1016/j.jembe.2010.03.002

70. Kenitz KM, Visser A, Ohman MD, Landry MR, Andersen KH. Community trait distribution across environmental gradients. Ecosystems. Ecosystems. 2019; 22, 968-980. https://doi.org/10.1007/s10021-018-0314-5

71. Taylor AG, Landry MR. Phytoplankton biomass and size structure across trophic gradients in the southern California Current and adjacent ocean ecosystems. Marine Ecology Progress Series. 2018; 592: 1-17. https://doi.org/10.3354/meps12526

72. Gorsky G, Ohman MD, Picheral M, Gasparini S, Stemmann L, Romagnan JB, et al. Digital zooplankton image analysis using the ZooScan integrated system. Journal of Plankton Research. 2010; 32: 285-303.

73. Ohman MD, Romagnan JB. Nonlinear effects of body size and optical attenuation on Diel Vertical Migration by zooplankton. Limnology and Oceanography. 2016; 61: 765-770. https://doi.org/10.1002/lno.10251

74. Lavaniegos BE, Ohman MD. Coherence of long-term variations of zooplankton in two sectors of the California Current System. Progress in Oceanography. 2007; 75: 42-69. https://doi.org/10.1016/j.pocean.2007.07.002

75. Venrick E. Floral patterns in the California Current: The coastal-offshore boundary zone, Journal of Marine Research. 2009; 67: 89-111.

76. Armstrong RA. Grazing limitation and nutrient limitation in marine ecosystems: Steady-state solutions of an ecosystem model with multiple food chains. Limnology and Oceanography. 1994; 39: 597-608. https://doi.org/10.4319/lo.1994.39.3.0597

77. Venrick E. Floral patterns in the California Current off southern California: 1990-1996, Journal of Marine Research. 2002; 60: 171-189.

78. Strom SL, Miller CB, and Frost BW. What sets lower limits to phytoplankton stocks in high-nitrate, low-chlorophyll regions of the open ocean? Marine Ecology Progress Series. 2000; 193: 19-31. https://doi.org/10.3354/meps193019

79. Olivares M, Saiz E, Calbet A. Ontogenetic changes in the feeding functional response of the marine copepod Paracartia grani. Marine Ecology Progress Series. 2019; 616: 25-35. https://doi.org/10.3354/meps12928

80. Jeong HJ, Yoo YD, Kim JS, Seong KA, Kang NS, Kim TH. Growth, feeding and ecological roles of the mixotrophic and heterotrophic dinoflagellates in marine planktonic food webs. Ocean Science Journal. 2010; 45: 65-91. https://doi.org/10.1007/s12601-010-0007-2

81. Paffenhofer GA, Stearns DE. Why is *Acartia tonsa* (Copepoda, Calanoida) restricted to nearshore environments? Marine Ecology Progress Series. 1988; 42: 33-38. https://doi.org/10.3354/meps042033

82. Lessard EJ, Murrell MC. Microzooplankton herbivory and phytoplankton growth in the northwestern Sargasso Sea. Aquatic Microbial Ecology. 1998; 16: 173-188. https://doi.org/10.3354/ame016173

83. Vidal J. Physioecology of zooplankton. I. Effects of phytoplankton concentration, temperature, and body size on the growth rate of *Calanus pacificus* and *Pseudocalanus* sp. Marine Biology. 1980; 56: 111-134.

84. Kiørboe T, Visser AP, Andersen KH. A trait-based approach to ocean ecology. ICES Journal of Marine Science. 2018; 75: 1849-1863. https://doi.org/10.1093/icesjms/fsy090

85. Steele JH, Henderson EW. The role of predation in plankton models. Journal of Plankton Research. 1992; 14: 157-172.
86. Ohman MD, Hirche H-J. Density-dependent mortality in an oceanic copepod population. Nature. 2001; 412: 638-641. https://doi.org/10.1038/35088068 PMID: 11493921

87. Goebel N, Edwards CA, Zehr JP, Follows MJ. An emergent community ecosystem model applied to the California Current System. Journal of Marine Systems. 2010; 83: 221-241. https://doi.org/10.1016/j.jmarsys.2010.05.002

88. Goebel N, Edwards CA, Zehr JP, Follows MJ, Morgan SG. Modeled phytoplankton diversity and productivity in the California Current System. Ecological Modelling. 2013; 264: 37-47. http://dx.doi.org/10.1016/j.ecolmodel.2012.11.008

89. Chenillat F, Franks PJS, Rivière P, Capet X, Grima N, Blanke B. Plankton dynamics in a cyclonic eddy in the Southern California Current System, J. Geophys. Res. Oceans. 2015; 120, 5566-5588. https://doi.org/10.1002/2015JC010826

90. Kahru M, Mitchell BG. Empirical chlorophyll algorithm and preliminary SeaWiFS validation for the California Current. International Journal of Remote Sensing. 1999; 20: 3423-3429.