An investigation of grazing behaviors that result in winter phytoplankton biomass accumulation

Recent observations have shown that phytoplankton biomass increases in the North Atlantic during winter, even when the mixed layer is deepening and light is limited. Current theories suggest that this is due to a release from grazing pressure. Here we demonstrate that the often-used grazing models that are linear at low phytoplankton concentration do not allow for a wintertime increase in phytoplankton biomass. However, certain mathematical formulations of grazing that are quadratic (or more generally non-linear) in phytoplankton concentration at low concentrations can reproduce the fall to spring transition 5 in phytoplankton, including wintertime biomass accumulation. We illustrate this point with a minimal model for the annual cycle of North Atlantic phytoplankton designed to simulate phytoplankton concentration as observed by BioGeoChemicalArgo (BGC-Argo) floats in the North Atlantic. This analysis provides a mathematical framework for assessing hypotheses of phytoplankton bloom formation.


Predator-prey decoupling
In this section we formulate a simple ecosystem model and examine different grazing functions to clarify the relationship between grazing rates and mixed layer depth during winter conditions (Figure 1).
Marine planktonic ecosystem dynamics can be coarsely represented as an interaction between three compartments: nutrients, phytoplankton, and zooplankton. These broad compartments integrate across all the chemical and biological diversity observed 60 in the ocean and are defined by their interactions with each other. In the simple formulation adopted here, the nutrients are consumed by phytoplankton, the zooplankton consume phytoplankton, and the plankton are converted back to nutrients when they die. The set of equations that describe these interactions for the concentrations of nutrients (n), phytoplankton (p), and zooplankton (z) as a function of the ocean depth ζ take the form: Dn Dt = −µ(n, t)e K d ζ p + d p p + (1 − a)g(p)z + d z z 2 + ∂ ∂ζ κ ∂n ∂ζ , Dp Dt = µ(n, t)e K d ζ p − g(p)z − d p p + ∂ ∂ζ κ ∂p ∂ζ , Dz Dt = ag(p)z − d z z 2 + ∂ ∂ζ κ ∂z ∂ζ .
(1) 65 The three compartments are modelled in terms of their carbon content and thus have the same units. The vertical coordinate, ζ, is zero at ocean surface and negative below. All compartments are mixed in the vertical by ocean turbulence at a rate set by the diffusivity κ. The phytoplankton specific growth rate depends on nutrients, according to the function µ(n), and decays exponentially with depth due to the absorption of light with depth with an attenuation coefficient K d . Phytoplankton mortality, −d p p, is linear in p. Zooplankton mortality, −d z z 2 , is quadratic in z to account for the drop in grazing by higher trophic levels 70 when food is scarce; this choice has the additional property of preventing extinction of zooplankton in winter. The grazing of phytoplankton by zooplankton is linear in z and proportional to p according to the grazing function g(p). Zooplankton are messy eaters and ingest only a fraction a < 1 of g(p)z.
We will focus on analyzing the model in equation 1 during the wintertime period through the bloom onset. During this period, we can make a few simplifying assumptions. First, we will assume that turbulence is strong enough to keep all compartments 75 well mixed in the vertical over a mixed layer of depth H. This assumption holds if the turbulence mixes all compartments throughout H on a timescale faster than any biological timescale (Taylor and Ferrari, 2011). Second, we will assume that winter growth is not nutrient limited (n n 0 ) and thus µ(n) saturates to a constant mixed layer-averaged value µ 0 (we will not make this assumption about a time and nutrient independent growth rate in section 3). For simplicity, we also assume that nutrients are exactly conserved in the mixed layer and model the nutrient concentration implicitly as n = N max − p − z, where 80 N max is the total carbon content which gets redistributed across the three compartments. Finally, we assume that the mixed layer is deep relative to the depth of light penetration (HK d 1). All these assumptions are appropriate in winter, the focus of our study, but they are less defensible in other seasons when turbulence is weak (Taylor and Ferrari, 2011).
We formulate a bulk mixed layer model by employing these simplifying assumptions and taking the vertical average of the equations in (1) over the mixed layer depth H(t), g0 and p0 are given in Table 1  where p and z are the constant mixed layer concentrations of phytoplankton and zooplankton, respectively. The term is the average light over the mixed layer, which is computed as the integral of the light level over the mixed layer depth divided by the mixed layer depth.
The term s + appears when taking the vertical average of the mixing term in equation (1). It represents the dilution of 90 phytoplankton and zooplankton that results from the turbulent entrainment of water without biomass across the mixed layer base and is given by, We can derive an equation for the standing stock of biomass in the mixed layer by taking a vertical integral of the equations in (1). Introducing P = Hp and Z = Hz to represent the total biomass of phytoplankton and zooplankton respectively we have 4 https://doi.org/10.5194/bg-2020-444 Preprint. Discussion started: 1 December 2020 c Author(s) 2020. CC BY 4.0 License.
In contrast to the average concentration, the total biomass does not change due to the physical effects of dilution. However, when the mixed layer shoals, biomass is lost below the mixed layer and the total biomass decreases at a rate given by s − In the following subsections we will analyze the phenology of phytoplankton for different choices of grazing functions (Figure 1). The linear (Holling type I) grazing function assumes that the plankton-specific grazing rate (units of per day) increases linearly with phytoplankton concentration, g H I (p) = g 0 p. The saturating functional responses are linear at low prey concentration and saturate at high prey concentrations. An example saturating response is the Holling type II functional response, . This functional response assumes that processing of food and searching for food are mutually exclusive 105 behaviors (Visser, 2007;Kiørboe et al., 2018). The parameter g 0 is a function of processing time and the parameter p 0 is a function of both search and processing time. This parsimonious theoretical basis and ability to fit the parameters from experimental data has made this functional response one of the most commonly used (Verity, 1991;Kiørboe et al., 2018). The Holling type III functional response has a reduction in grazing at low prey concentration. One formulation is a sigmoidal function, g H III (p) = g 0 p 2 p 2 0 +p 2 which is quadratic in p for low p. This type III functional response can be derived as a generalization 110 of the type II response where the search time is a linear function of prey concentration. Effectively, there is a prey refuge at low concentration because it is more difficult for predators to find each prey item. There are other possible mechanisms for a type III functional response, including a threshold response by predators (Mullin et al., 1975;Ohman, 1984) and prey switching (Vallina et al., 2014). To compare the functional responses, we formulate the zooplankton specific grazing rate as a function of phytoplankton concentration 115 g(p) = g 0 (p/p 0 ) k−1 1 + (p/p 0 ) k−1 .
The exponent k determines the degree of non-linearity of the functional response. The Holling type II functional response is k = 2 and the Holling type III response is k = 3. The parameter p 0 is a half saturation constant. When p = p 0 , the grazing is at half of the maximum rate g(p 0 ) = g0 2 .
2.1 Grazing linear in phytoplankton concentration for constant zooplankton concentration: Sverdrup's model 120 Phytoplankton are known to respond faster than zooplankton to environmental changes (Fileman and Leakey, 2005). Sverdrup (1953) assumed that such an assumption applied to the rapid onset of the spring bloom and proposed to model the phytoplankton growth rate according to equation (2), but setting g(p)z = g0 p0 z 0 p with z 0 the constant zooplankton concentration before the bloom onset, 125 Sverdrup focused on the time at the end of winter when the mixed layer starts shoaling in response to spring atmospheric conditions and thus could ignore the entrainment, i.e. s + = 0. Under these assumptions the mixed layer depth H(t) is the only time dependent parameter which can determine whether the phytoplankton concentration is exponentially decaying (winter conditions) or exponentially increasing (spring bloom onset). This gave rise to the widely applied 'critical depth hypothesis' which states that phytoplankton accumulation starts when the mixed layer shoals beyond a critical depth, While the critical depth hypothesis has become the most widely accepted framework to interpret the onset of spring blooms-but there are growing objections Behrenfeld (2010)-it is not very useful to make quantitative predictions. The criterion requires knowledge of the grazing rate at the end of winter before bloom onset, which is very difficult to measure. Sometimes this obstacle is overcome by assuming that g 0 z 0 d p , in which case the critical depth dependence on grazing can be ignored.

135
However the assumption is likely inappropriate for most blooms. For example, assuming a typical attenuation coefficient of K d = 0.05 m −1 in the winter North Atlantic (Organelli et al., 2017), where bloom onset is often observed at a critical depth of around 200 m (as reported in Siegel et al. (2002)), the ratio of growth to mortality rate, dp+g0z0/p0 µ0 , is predicted to be close to be 0.1. Mortality timescales of phytoplankton are believed to be much longer than ten times their division rates implying that grazing, not mortality, dominates phytoplankton losses at bloom onset. A theory of blooms must therefore 140 include a predictions of the zooplankton concentrations and their grazing rates at the end of winter, if it is to make falsifiable predictions. Additionally, on seasonal timescales there is substantial variation in zooplankton concentrations so a theory that includes variable phytoplankton and zooplankton concentrationss is necessary. The goal of the next two sections is to present two models of grazing with a focus on wintertime conditions.

Grazing linear in phytoplankton concentration
Consider first the saturating (type II) grazing function. In winter, prey concentrations are very low and this function is approx- Figure 2b). During the wintertime, as the mixed layer deepens, water from below the mixed layer is entrained, decreasing the concentration of the phytoplankton and zooplankton (s + > 0) but not their standing stock Assuming the natural mortality of phytoplankton is negligibly small, the growth and grazing terms in the P and Z equations have the same dependence on mixed layer depth H and thus any increase in H does not reduce grazing any more than it reduces the growth of phytoplankton. Consider for example a population in equilibrium, i.e. dP dt = dZ dt = 0. The equilibrium populations are Z * ≈ µ0p0 are independent of H, and thus an equilibrium population will remain in equilibrium even as the mixed layer deepens (Figure 2a). If phytoplankton biomass decreases at some point in winter then 155 subsequent changes in mixed layer depth cannot trigger any biomass accumulation as long as the biological parameters µ 0 , a, g 0 , d p , d z , and K d remain constant.
It could be rebutted that winter accumulation is possible if zooplankton mortality is represented as a linear, rather than quadratic loss term. In that case, as the mixed layer deepens, zooplankton biomass loss rates would not decrease as quickly as in the grazing rate as the mixed layer increases. This occurs at low values of phytoplankton biomass and deep mixed layers, as required for the disturbance-recovery hypothesis. The grazing rates with Holling type I functional response is independent of mixed layer depth. At high biomass there is also a decrease in grazing rate at both deeper and shallower mixed layers with a type III functional response and at shallow mixed layers with a type II functional response.
the rate of zooplankton grazing on phytoplankton, eventually reaching a crossing over point at which there would be a net loss 160 of zooplankton biomass and consequently an increase in phytoplankton biomass. This is the case of Lotka-Volterra predatoryprey dynamics in a variable environment (Yorke and Anderson Jr, 1973;Dubois, 1975). However, this model is problematic because a linear zooplankton mortality at low concentrations is only defensible in the absence of grazing by higher trophic levels. Such grazing is what is implicitly modelled with a quadratic mortality term such as the one used in equation (9).

Grazing quadratic in phytoplankton concentration
The situation is different if we prescribe a phytoplankton grazing function with a dependence on p that is stronger than linear.
The Holling type III functional response is a popular choice and can be written as g H III (p) = g 0 p 2 p 2 +p 2 0 , which can be approximated as g H III (p) ≈ g0 p 2 0 p 2 at low prey concentration. With this functional response, the rate of change of biomass is given by, 170 In this case, the grazing rate decreases faster than the phytoplankton growth rate as the mixed layer deepens due to the additional 1 H factor in the grazing term ( Figure 2c). This opens the possibility of a net increase in phytoplankton biomass due to deepening of the mixed layer, consistent with the scenario invoked in the disturbance-recovery hypothesis (Behrenfeld, 2010;Behrenfeld and Boss, 2014). This is the key result of this paper. In what follows, we will use observations to explore the implications of this insight beyond the low phytoplankton limit.
We aim to demonstrate that when implemented in a full NPZ model, a grazing function with a quadratic (or higher) dependence on phytoplankton concentration at low p is sufficient to reproduce both wintertime biomass accumulation and a spring bloom. In order to model the full annual cycle, we utilize a more realistic phytoplankton growth rate that depends on nutrient concentration and has a yearly cycle. We replace the growth term µ0 This growth rate has temporal dependence through the mixed layer depth and through the surface irradiance. It also depends on nutrient concentration through the function n/(n + n 0 ) which varies throughout the year and is close to one in winter when p and z are small. Using this formulation we now test the impact of the grazing function on the yearly evolution of biomass and compare with in-situ observations.

185
Float measurements of the phytoplankton annual cycle in the North Atlantic We calibrate the NPZ model using the averaged annual cycle of phytoplankton biomass as observed by BGC-Argo floats in the high-latitude North Atlantic (Mignot et al., 2018). Our model ignores the effect of lateral heterogeneity or restratification on phytoplankton dynamics. In order to relate the model results to empirical data, we followed Mignot et al. We estimated phytoplankton concentration p(t) from backscatter measurements. The mixed layer depth H(t) is defined as the depth at which the potential density increases by 0.03 kg m −3 from the potential density at 10m. As in Mignot et al. (2018) 195 the net phytoplankton population growth rate was then calculated using the observed phytoplankton concentration and mixed layer depth as In contrast to Mignot et al. (2018), the accumulation rate was computed over the mixed layer rather than the productive layer.
In order to account for interannual and regional variability in bloom timing, we rescaled the time axis of each individual float

Model parameters
The NPZ model equations (2) are solved replacing µ 0 with µ(n, t) as given in equation (11) and using the yearly timeseries of H(t) estimated from the average from all float measurements. The total nitrogen pool is assumed to be constant at N max = 30 mg C m −3 which implies that phytoplankton and zooplankton are immediately remineralized once they die so 210 that n + p + z = N max at all times and one needs equations only for p and z, while the nutrient concentration can be inferred from n = N max − p − z. Some parameters are prescribed based on reasonably well established values found in the literature: (Eppley, 1972;Geider et al., 1998;Bissinger et al., 2008), a=0.5 (Landry et al., 1984;Moore et al., 2001), et al., 2001), and K d =0.05 m −1 . However other parameters relating to grazing and zooplankton and phytoplankton mortality are much more uncertain (see Table 1). These parameters and the initial p and z concentrations are 215 therefore calibrated by fitting the model accumulation rate and concentrations to observations over the full annual cycle using a trust-region-reflective least-squares algorithm (Coleman and Li, 1996). The accumulation rates are smoothed before fitting with a five-point Savitzky-Golay filter. The best fit parameters values are given in Table 1. The temporal rescaling used to average the timeseries creates a spurious peak in net population growth rate at the beginning of winter. Throughout the winter there is variability in accumulation rates among individual timeseries, including some negative

Comparison of model and observations
Using either Holling type II or III grazing functions, the model with the best fit parameters generates a spring bloom with a rapid increase in phytoplankton concentration and biomass that coincides with the spring shoaling of the mixed layer ( Figure   230 3b). However, the Holling type III model results in net positive phytoplankton population growth through the winter, while the Holling type II model does not (Figure 3a). The commonly used grazing functions of Holling type I and II do not satisfy the requirement of non-linear dependence of grazing on phytoplankton at low phytoplankton concentrations, and thus cannot capture the observed wintertime biomass accumulation, while the Holling type III functional response has the appropriate nonlinear dependence (Holling, 1959).

235
During the winter, the phytoplankton concentration is larger when using the type III grazing function than when using the type II grazing function (Figure 3b). Despite this, the winter grazing rate is lower with the type III grazing function (Figure 3a).
In order to compensate for the larger winter grazing with the type II function, the parameter fitting procedure infers a much lower linear phytoplankton mortality d p for that case (Table 1).
During the summertime, the mixed layer depth is fairly constant and phytoplankton and zooplankton populations are close 240 to equilibrium. This model does not include export from the mixed layer through sinking or migrating particles. Instead, any carbon export from the mixed layer only occurs when the mixed layer is shoaling due to biomass being left in the stratified layer below the new mixed layer.
The modeled relationship between phytoplankton and zooplankton shows notable differences between the two grazing functions. This is best illustrated by plotting the temporal evolution of the two communities in a z − p plane as shown in Figure 4.

245
At the end of winter, zooplankton are at slightly higher concentration with type III than type II grazing, because they have fared better throughout the winter by feeding on a larger phytoplankton population. Zooplankton respond slowly to the explosive spring phytoplankton bloom with the type II grazing, resulting in higher phytoplankton growth rates and a lower zooplankton concentration. By contrast, with the type III grazing, zooplankton are strongly coupled to phytoplankton and start grazing as soon as the bloom gets going reducing its amplitude. Importantly the rate of increase of phytoplankton concentration during 250 the spring bloom is slower than exponential (Mignot et al., 2018), consistent with the prediction of the disturbance-recovery hypothesis (Behrenfeld and Boss, 2014). With both grazing functions, the spring bloom populations are out of equilibrium with phytoplankton concentrations being higher and zooplankton concentrations being lower than at equilibrium.
The simple n − p − z model used here is an imperfect representation of the observations. For example, the model only includes one phytoplankton type and one zooplankton type, which precludes both the succession of different phytoplankton 255 types during the spring and summer and the presence of a microbial loop that could reduce the flow of carbon up the food chain (Azam et al., 1983). The bulk zero-dimensional model assumes that phytoplankton and zooplankton concentrations are uniform in the mixed layer and zero below, a defensible approximation for winter conditions-the focus of this study-when the mixed layer is deep and turbulent mixing is strong, but not in other seasons when turbulence is weak (Taylor and Ferrari, 2011). The deficiencies of the bulk model are evident at the spring bloom onset, which is slightly delayed in the model relative 260 to the observations, occurring once the mixed layer has shoaled rather during mixed layer shoaling. In observations, blooms start as soon as turbulent mixing subsides because phytoplankton is no longer mixed away from the surface, while there is a lag of days to weeks before the mixed layer restratifies and shoals (Taylor and Ferrari, 2011). The bulk model is also problematic in summer when the mixed layer is shallower than the euphotic layer and some of the productivity takes place below the mixed layer base where the model assumes p = z = 0. Despite these deficiencies, bulk mixed layer models have been shown 265 to qualitatively reproduce the full annual cycle of plankton dynamics in other regions (c.f. Evans and Parslow (1985)) and are especially appropriate for our work which focuses on phytoplankton growth in winter.

Discussion
Our work suggests that the winter accumulation of biomass recently documented from float observations in the North Atlantic (Behrenfeld, 2010;Mignot et al., 2018), while much weaker than that the spring and summer accumulation (Lutz et al., 270 2007;Uitz et al., 2010), reveals otherwise hard to document top-down controls on phytoplankton populations. By studying winter time phytoplankton population dynamics, when growth conditions are less than optimal, we have been able to make inferences about the rate of zooplankton grazing. We demosntrated that there must be non-linearity in the interaction between the zooplankton and the phytoplankton at low phytoplankton concentrations in order to release the phytoplankton from grazing pressure. A quadratic grazer response function at low phytoplankton biomass is sufficient for phytoplankton biomass accumu-275 lation, although higher order nonlinearities would also reproduce the observed dynamics. Relatively little is known about phytoplankton loss through grazing (Dolan and McKeon, 2005) in comparison to the other factors that control the dynamics of phytoplankton populations like macro-and micronutrients, light availability, and temperature (Eppley, 1972). Our work suggests that winter conditions may offer a unique opportunity to study phytoplankton grazing in the field. In the wintertime, cell division rates decrease because of light limitation due to both deepening of the mixed layer 280 and decrease in sea surface light. In order for phytoplankton accumulation rates to be positive while cell division rates are declining, the phytoplankton loss rates must decrease faster than division rates. While grazing is not the only concentrationdependent process, it is an interesting and compelling example of one process that could lead to biomass accumulation during the wintertime and it is the fundamental tenet of the "disturbance-recovery hypothesis" (Behrenfeld, 2010).
The wintertime growth is not only important to sustain phytoplankton populations in winter, but it is believed to play a 285 crucial role for the development of the subsequent spring bloom. We showed that a non-linear coupling between zooplankton and phytoplankton results in larger wintertime populations of both zooplankton and phytoplankton at the end of winter than would occur with a linear coupling. Furthermore a non-linear coupling results in a faster acceleration in zooplankton grazing once phytoplankton concentrations increase during a bloom. The combination of more abundant wintertime populations and stronger/more rapid coupling between phytoplankton and zooplankton populations curb explosive phytoplankton growth.

290
There are other possible explanations for wintertime biomass accumulation beyond the dilution of phytoplankton. The biological functions encapsulated in the NPZ model parameters may vary over time. For example, the zooplankton assimilation rate, a, could depend on the nutrient content of the prey introducing an alternative nonlinear effect (Landry et al., 1984;Irigoien et al., 2005). In our model, the only time-dependent terms are light availability, which has little influence on wintertime biomass accumulation, and mixed layer depth, which drives the wintertime biomass accumulation. Finally, there is evidence 295 that wintertime growth can be triggered by mixed layer instabilities that occasionally restratify the mixed layer during the winter and thus increase the light available for phytoplankton (Karimpour et al., 2018). However this cannot be the unique explanation, because float observations presented in (Mignot et al., 2018) and reviewed here show many examples of wintertime accumulation where these mixed layer dynamics did not seem to apply.
The sensitive dependence of phytoplankton phenology on the rate of grazing by higher trophic levels at low concentrations 300 provides a powerful quantitative framework in which to evaluate theories of plankton phenology. Observations of wintertime phytoplankton biomass accumulation have been interpreted as evidence of a release from grazing pressure in deep mixed layers, but little attention has been given to the key role played by the choice of grazing functions in these theories. Some studies have used a Holling type III grazing function (Behrenfeld and Boss, 2014;Yang et al., 2020) while others relied on a prey switching formulation, where the zooplankton preferentially consumes the most common type of phytoplankton (Llort 305 et al., 2015). The observational evidence for a lower bound on phytoplankton concentrations (Lessard and Murrell, 1998) ought to be studied within the framework presented here. The interpretation of the response of phytoplankton to sudden environmental perturbations on subseasonal timescales, such as storms (Behrenfeld and Boss, 2018), will also require a careful assessment of the grazing functions which control how fast zooplankton grazing responds to increases in phytoplankton concentrations.
Last, but not least, this framework ought to be applied to test the predictions of different theories of bloom onset (Verity et al.,310 1993; Morison et al., 2020;Mojica et al., 2020). While our work pointed out the key role played by the choice of grazing functions in such theories and models, it is important to point out that state of the art Earth System Models often use multispecies ecosystem models (Laufkötter et al., 2015). Multispecies models do not necessarily result in the same dynamics as the single-species functional responses used here (Gentleman et al., 2003), but our result that phytoplankton phenology is very sensitive to the degree of non-linearity in growth and mortality functions is very likely to hold for more complex ecosystem 315 models as well.
It is worth commenting on the ecological underpinnings for the different models of grazing. The grazing functions used in our model represent the coupling between all species of each trophic level of phytoplankton and zooplankton; the phytoplankton class includes all autotrophs, while the zooplankton class includes all grazers that consume phytoplankton. A nonlinear decrease in grazing rates at low prey concentration has been observed in the lab studies of aquatic vertebrates and invertebrates and in 320 theoretical studies (Real, 1977;Barrios-O'Neill et al., 2016), but it is not clear how these results translate to full populations.
Similarly, copepods go into diapause in the wintertime (Baumgartner and Tarrant, 2017), and this effectively reduces grazing pressure during winter; however, microzooplankton account for the majority of phytoplankton mortality in the ocean (Landry and Calbet, 2004). Other mechanisms described by the Holling type III functional response include prey switching (Vallina et al., 2014), predator learning (Holling, 1966), and prey refuges (Taylor, 2013). While there are few structural prey refuges in 325 an oceanic mixed layer, a patchy environment can also provide a type of prey refuge. A Holling type III functional response can arise due to non-random grazing behavior when population dynamics are integrated over a patchy environment (Nachman, 2006;Morozov, 2010).

Conclusions
A reduction in the grazing rate at low phytoplankton concentration has been proposed as the mechanism to explain the emerging 330 observation that biomass often increases, albeit weakly, during the wintertime when mixed layers deepen (Behrenfeld, 2010).
It has also been pointed out that the critical depth hypothesis fails to capture wintertime growth because it implicitly assumes that loss rates are constant either because they are dominated by constant respiration or by grazing by a constant zooplankton population. Previous modeling results have not acknowledged that a reduction in grazing pressure through dilution of plankton populations in deep mixed layers requires a grazing function that decreases faster than linearly in phytoplankton concentration  (2005). This coupling has been interpreted as evidence that zooplankton grazing pressure can respond very rapidly to any 340 changes in phytoplankton concentrations (Behrenfeld and Boss, 2018). To the extent that these interpretations are correct, our study suggests that observations can therefore be used to infer the functional form of grazing functions, an aspect of plankton ecology that is otherwise very to quantify.
Observational validation of the functional forms of grazing functions is key to build confidence in predictions based on biogeochemical models. Different models can be tuned to provide reasonable estimates of the annual cycle of phytoplankton 345 biomass, like our NPZ model with a more linear grazing function. However, in order to make predictions that are robust to changing conditions, it is important that models have the correct functional dependencies. Tuning of model parameters is no guarantee of model performance in an evolving environment that has not been observed yet. Climate change may reshape North Atlantic phytoplankton populations and primary production (Balaguru et al., 2018) due to increasing surface temperature, shoaling mixed layer depths, and increasing upper ocean stratification (Edwards and Richardson, 2004). Predicting and 350 quantifying the impact of these changes requires robust model formulations, not models tuned to present climate conditions.
Code and data availability. All code and data are available at https://github.com/mara-freilich/grazing_functions_bg