Reliable underwater acoustic propagation is challenging due to complex ocean dynamics such as internal-waves and to the uncertain larger-scale ocean physics, acoustics, bathymetry, and seabed fields. For accurate acoustic propagation, capturing the important environmental uncertainties and variabilities and predicting the probability distributions of the acoustic pressure field is then what matters. Prior works towards addressing this goal include (i) wave propagation in random media techniques such as perturbation methods, path integral theory, and coupled-mode transport theory, and (ii) probabilistic modeling techniques such as Monte Carlo sampling and Polynomial Chaos expansions. Recently, we developed a novel technique called the Dynamically Orthogonal Parabolic Equations (DO-ParEq) which represent the sound speed, density, bathymetry, and acoustic pressure fields using optimal dynamic Karhunen-Loeve decompositions. The DO-ParEq are range-evolving partial and stochastic differential equations preserving acoustic nonlinearities and non-Gaussian properties. In this presentation, we showcase the theoretical and computational advantages of the DO-ParEq framework compared to the state-of-the-art techniques in the Pekeris waveguide and wedge benchmark problems, in addition to a realistic ocean example in the New York Bight region.