More accurate methods are needed to represent biogeochemistry in ocean models with coarse spatial resolution, in order to assess the response of marine ecosystems to global change. We use eddy-resolving simulations to test methods of upscaling biogeochemistry from 1 km to the 100 km scale of global model grid cells. The neglect of subgrid-scale variability results in serious errors which are not robustly corrected by retuning parameters in the model dynamics. Moment closure schemes provide accurate upscaling for modest computational investment, with broadly similar results obtained by second moment and conditional moment closure schemes. However, the conditional scheme gives clear improvement when variability is imposed on maximum uptake rates under Michaelis–Menten nutrient limitation, as this may invalidate second-order expansions of the mean field dynamics.