A 4D, multivariate probability distribution is derived of plankton ecosystem variables on Georges Bank given the data from the GLOBEC program. Covariance models for chlorophyll, nutrient, zooplankton, temperature and surface irradiance fields are fitted by maximum likelihood estimation, allowing for anisotropy, non-stationarity, cross-correlations and distributional uncertainty within the power-normal family. Model selection criteria are applied and the fitted models are tested and bias-corrected using parametric bootstrap methods. The probability model not only serves as a means of interpolating the data, but also allows us to simulate plausible configurations of the fields (conditional simulation). Hence we can study the effects of errors in forcings and initial/boundary conditions in spatial and non-spatial dynamical plankton models, after fitting the biological model parameters. The simulations also enable us to test plankton model formulations against the data whilst allowing for these residual sources of error.