Predictive dynamical models for marine ecosystems are used for a variety of needs. Due to the sparse measurements and limited understanding of the myriad of ocean processes, there is however significant uncertainty. There is model uncertainty in the parameter values, functional forms with diverse parameterizations, and level of complexity needed, and thus in the state variable fields. We develop a Bayesian model learning methodology that allows interpolation in the space of candidate dynamical models and discovery of new models from noisy, sparse, and indirect observations, all while estimating state variable fields and parameter values, as well as the joint probability distributions of all learned quantities. We address the challenges of high-dimensional and multidisciplinary dynamics governed by partial differential equations (PDEs) by using state augmentation and the computationally efficient Gaussian Mixture Model – Dynamically Orthogonal filter. Our innovations include stochastic formulation parameters and stochastic complexity parameters to unify candidate models into a single general model as well as stochastic expansion parameters within piecewise function approximations to generate dense candidate model spaces. These innovations allow handling many compatible and embedded candidate models, possibly none of which are accurate, and learning elusive unknown functional forms that augment these models. Our new Bayesian methodology is generalizable and interpretable. It seamlessly and rigorously discriminates among existing models, but also extrapolates out of the space of models to discover new ones. We perform a series of twin experiments based on flows past a ridge coupled with three-to-five component ecosystem models, including flows with chaotic advection. We quantify the learning skill, and evaluate convergence and the sensitivity to hyper-parameters. Our PDE framework successfully discriminates among functional forms and model complexities, and learns in the absence of prior knowledge by searching in dense function spaces. The probabilities of known, uncertain, and unknown model formulations, and of biogeochemical-physical fields and parameters, are updated jointly using Bayes’ law. Non-Gaussian statistics, ambiguity, and biases are captured. The parameter values and the model formulations that best explain the noisy, sparse, and indirect data are identified. When observations are sufficiently informative, model complexity and model functions are discovered.