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

•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.•The closure modeling approach is used to account for spatial variability and temporal fluctuations (higher central moments) of the phytoplankton distribution.•NP, NPZ and NPZD plankton models are developed 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 the NPZ model, including one predator, can reproduce consistently the full range of observed variability, as quantified by the observable CVp, even under different assumptions about the functional dependence of mortality.•For consistently making testable predictions about the response of plankton ecosystems under micro-scale variability, models need to include at minimum one class of predators (Z, herein). 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