A new methodology for rigorous Bayesian learning of high-dimensional stochastic dynamical models is developed. The methodology performs parallelized computation of marginal likelihoods for multiple candidate models, integrating over all state variable and parameter values, and enabling a principled Bayesian update of model distributions. This is accomplished by leveraging the dynamically orthogonal (DO) evolution equations for uncertainty prediction in a dynamic stochastic subspace and the Gaussian Mixture Model-DO filter for inference of nonlinear state variables and parameters, using reduced-dimension state augmentation to accommodate models featuring uncertain parameters. Overall, the joint Bayesian inference of the state, model equations, geometry, boundary conditions, and initial conditions is performed. Results are exemplified using two high-dimensional, nonlinear simulated fluid and ocean systems. For the first, limited measurements of fluid flow downstream of an obstacle are used to perform joint inference of the obstacle’s shape, the Reynolds number, and the O(10⁵) fluid velocity state variables. For the second, limited measurements of the concentration of a microorganism advected by an uncertain flow are used to perform joint inference of the microorganism’s reaction equation and the O(10⁵) microorganism concentration and ocean velocity state variables. When the observations are sufficiently informative about the learning objectives, we find that our posterior model probabilities correctly identify either the true model or the most plausible models, even in cases where a human would be challenged to do the same.