Understanding plankton ecosystem dynamics under realistic micro-scale variability requires modeling at least three trophic levels

Plankton ecosystems are complex, multi-trophic networks of biotic and abiotic interactions among physical and chemical components. Still, nutrient-phytoplankton (NP) interactions are in many cases assumed to be representative of higher trophic-level interactions in plankton ecosystems. Here, we investigate the degree to which NP interactions capture the overall dynamics of multi-trophic plankton ecosystems when accounting for realistic levels of micro-scale variability. Plankton models are typically developed based on the mean-field approach, which considers only first central moments (i.e., spatio-temporal means). Such conventional plankton models may be appropriate for meso‑ or larger-scales, but inappropriate for the highly intermittent spatial fluctuations of phytoplankton that are ubiquitous at the micro (mm) scale. Using Reynold's decomposition, the closure approach accounts for spatial variability and temporal fluctuations (higher central moments) of the phytoplankton distributions. We apply closure models of various combinations of Nutrient (N), Phytoplankton (P), and Zooplankton (Z), each with linear, hyperbolic, sigmoidal and quadratic phytoplankton mortality forms to test i) the previously advanced hypothesis ‘fluctuations enhance higher trophic level biomass in plankton models’ and ii) the suitability of nutrient-phytoplankton interactions as a proxy for higher trophic ecological processes. We find that in NP models with linear, hyperbolic and sigmoidal phytoplankton mortality (P-mortality), phytoplankton biomass (P-biomass) increases with variability, but in the NP model with quadratic P-mortality, P-biomass decreases with variability. The results are robust in the sense that similar qualitative behavior is obtained for a wide range of biologically feasible parameter sets. We investigate the characteristics of the measurable quantity coefficient of variation of phytoplanktonCVp for models with different numbers of compartments. The NP models with density-independent form (linear form) and density-dependent forms (hyperbolic, sigmoidal and quadratic) for P-mortality term behave quite opposite to each other: NP models with linear P-mortality produce coefficient of variation CVp > 1 for all stable solutions while NP models with hyperbolic, sigmoidal and quadratic P-mortality forms produce CVp < 1 for all stable solutions. Depending on the mortality rate assumed, the NP model produces qualitatively different results at both micro- and macro-scales, and therefore is unable to consistently capture the true underlying dynamics of plankton ecosystems. By contrast, NPZ models behaved consistently and supported the previously advanced hypotheses, regardless of the functional form assumed for P-mortality.