The last two decades have seen the development of a few high-order nonhydrostatic (NHS) ocean models, but there is still a need to better understand the associated numerical properties. A key improvement in NHS models in comparison to hydrostatic models is their dispersive nature, which allows for accurate resolution of gravity waves. However, it has been shown (Vitousek 2011) that low-order finite difference and finite volume methods can suffer from large numerical dispersion, which hinders the resolution of such waves. Therefore, high-order methods pose an attractive alternative, and in this study, we explore the use of high-order hybridizable discontinuous Galerkin-based (HDG) schemes (Ueckermann 2016) to resolve NHS gravity waves. We quantify and compare the numerical dispersion and computational cost of these schemes to their low-order counterparts using idealized internal wave and bottom gravity current test cases. Additionally, the stability of NHS models in the presence of fast-moving free-surface gravity waves is crucial to their efficacy. To this end, we explore the stability and convergence of high-order HDG methods in the context of the NHS primitive ocean equations solved on skewed computational domains. Finally, using our high-order HDG NHS solver, we illustrate results from process studies of gravity-driven wave dynamics including (i) 3D internal wave propagation over complex bathymetry, (ii) tidally-forced oscillatory flow over seamounts and (iii) bottom gravity currents.