Reliable acoustic predictions remain challenging due to the sparse and heterogeneous data, as well as to the complex ocean physics, sea surface and seabed processes, multiscale interactions, and large dimensions.These complexities lead to several sources of uncertainty. Predicting the full probability distributions of the ocean-acoustic-seabed fields then allows robust informative modeling, inference, and decision-making. In this work, we integrate our acoustic stochastic Dynamically-Orthogonal Parabolic Equations (DO-PEs) and Gaussian Mixture Model-DO (GMM-DO) frameworks with the MSEAS primitive equation ocean modeling system to enable unprecedented probabilistic forecasting and learning of ocean physics and acoustic pressure and transmission loss (TL) fields, accounting for uncertainties in the ocean, acoustics, bathymetry, and seabed fields. We demonstrate the use of this system for low to mid-frequency propagation with real ocean data assimilation in three regions. The first sea experiment takes place in the western Mediterranean Sea where we showcase the system’s performance in predicting ocean and acoustic probability densities, and assimilating sparse TL and sound speed data for joint ocean physics-acoustics-source depth inversion in deep ocean conditions with steep ridges. In the second application, we simulate stochastic acoustic propagation in Massachusetts Bay around Stellwagen Bank and use our GMM-DO Bayesian inference system to assimilate TL data for acoustic and source depth inversion in shallow dynamics with strong internal waves. Finally, in the third experiment in the New York Bight, we employ our system as a novel probabilistic approach for broadband acoustic modeling and inversion. Overall, our results mark significant progress toward end-to-end ocean acoustic systems for new ocean exploration and management, risk analysis, and advanced operations.