Accurate sea ice models are essential to predict the complex evolution of rapidly changing sea ice conditions and study impacts on climate, wildlife, and navigation. However, numerical models for sea ice contain various uncertainties associated with initial conditions and forcing (wind, ocean), as well as with parameter values and parameterizations, functional forms of the constitutive relations, and state variables such as sea ice thickness and concentration, all of which limit predictive capabilities. In this work, we first develop new stochastic partial differential equation (PDE)-based Sea Ice Dynamically Orthogonal equations and schemes for efficient uncertainty propagation and probabilistic predictions. These equations and schemes preserve nonlinearities in the underlying spatiotemporal dynamics and evolve the non-Gaussianity of the statistics with a lower computational cost than Monte Carlo methods commonly used in sea ice data assimilation and sensitivity analysis. We then use the Gaussian Mixture Model (GMM)-DO filter for sea ice Bayesian nonlinear data assimilation and learning. Assimilating noisy and sparse measurements, we provide posterior probability distributions for not only the sea ice velocities, thickness, and concentration, but also for the external forcing, parameters, and even functional forms of the sea ice model. The equations and schemes are evaluated using stochastic test cases, in which we showcase the ability to evolve non-Gaussian statistics and capture complex nonlinear dynamics efficiently. We demonstrate the stochastic convergence of the probabilistic predictions to the stochastic subspace size and coefficient samples. Finally, we highlight the principled joint nonlinear inference and learning of the sea ice state and dynamics.