Retrospective inference through Bayesian smoothing is indispensable in geophysics, with crucial applications in ocean and numerical weather estimation, climate dynamics, and Earth system modeling. However, dealing with the high-dimensionality and nonlinearity of geophysical processes remains a major challenge in the development of Bayesian smoothers. Addressing this issue, a novel subspace smoothing methodology for high-dimensional stochastic fields governed by general nonlinear dynamics is obtained. Building on recent Bayesian filters and classic Kalman smoothers, the fundamental equations and forward–backward algorithms of new Gaussian Mixture Model (GMM) smoothers are derived, for both the full state space and dynamic subspace. For the latter, the stochastic Dynamically Orthogonal (DO) field equations and their time-evolving stochastic subspace are employed to predict the prior subspace probabilities. Bayesian inference, both forward and backward in time, is then analytically carried out in the dominant stochastic subspace, after fitting semiparametric GMMs to joint subspace realizations. The theoretical properties, varied forms, and computational costs of the new GMM smoother equations are presented and discussed.