Ocean plankton models are useful tools for understanding and predicting the behaviour of planktonic ecosystems. However, when the regions represented by the model grid cells are not well mixed, the population dynamics of grid cell averages may differ from those of smaller scales (such as the laboratory scale). Here, the 'mean field approximation' fails due to 'biological Reynolds fluxes' arising from nonlinearity in the fine-scale biological interactions and unresolved spatial variability. We investigate the domain-scale behaviour of two-component, 2D reaction-diffusion plankton models producing transient dynamics, with spatial variability resulting only from the initial conditions. Failure of the mean field approximation can be quite significant for sub grid-scale mixing rates applicable to practical ocean models. To improve the approximation of domain-scale dynamics, we investigate implicit spatial resolution methods such as spatial moment closure. For weak and moderate strengths of biological nonlinearity, spatial moment closure models generally yield significant improvements on the mean field approximation, especially at low mixing rates. However, they are less accurate given weaker transience and stronger nonlinearity. In the latter case, an alternative 'two-spike' approximation is accurate at low mixing rates. We argue that, after suitable extension, these methods may be useful for understanding and skillfully predicting the large-scale behaviour of marine ecosystems.