A new methodology for Bayesian inference of stochastic dynamical models is developed. The methodology leverages the dynamically orthogonal (DO) evolution equations for reduced-dimension uncertainty evolution and the Gaussian mixture model DO filtering algorithm for nonlinear reduced-dimension state variable inference to perform parallelized computation of marginal likelihoods for multiple candidate models, enabling efficient Bayesian update of model distributions. The methodology also employs reduced-dimension state augmentation to accommodate models featuring uncertain parameters. The methodology is applied successfully to two high-dimensional, nonlinear simulated fluid and ocean systems. Successful joint inference of an uncertain spatial geometry, one uncertain model parameter, and 0(105) uncertain state variables is achieved for the first. Successful joint inference of an uncertain stochastic dynamical equation and 0(105) uncertain state variables is achieved for the second. Extensions to adaptive modeling and adaptive sampling are discussed.