The simulation of nonlinear internal waves (NLIWs) is a challenging task for several reasons. The inherent long width of NLIWs and their long propagation distances necessitate high-aspect-ratio, anisotropic and usually deformed (due to complex bathymetry) computational domains. Furthermore, the broad range of scales for turbulence-resolving simulations pose a great challenge in environmental flow solvers. This dissertation discusses the mathematical formulation and numerical implementation in a high-performance-computing context of a hybrid Spectral Element Method and Fourier Galerkin (SEM/FG) flow solver specifically tailored to address these difficulties. In the first part of this work, an emphasis is given on the most computationally challenging implicit stage of the flow solver, the numerical solution of the domain-decomposed pressure Poisson equation (PPE). The resulting Schur complement system of equations of the PPE is iteratively solved using a deflated block-Jacobi preconditioned conjugate gradient solver which ensures a fast convergence rate independent of the domain's aspect-ratio, the number of elements in the along-wave propagating direction and the polynomial order per spectral element. Other implementation details regarding the solver's development are also reported including benchmarks of increasing complexity. In the next part of this dissertation, the approximate enforcement of a free-slip boundary condition (BC) in the context of the propagation of a NLIW over a non-zero-curvature bottom boundary is analyzed. It is shown that when an approximate free-slip BC is imposed, the velocity components can be treated as scalar quantities allowing the use of the same computational machinery as for the calculation of the density and the pressure. The effective numerical drag produced by the approximate free-slip BC is quantified thereby enabling the assessment of the accuracy of the tested approximations. Finally, in the last chapter of this dissertation, the SEM/FG flow solver development culminates with the simulation of a three-dimensional shoaling of a mode-one internal solitary wave (ISW) of depression. In these massively parallel, high-resolution and high-accuracy simulations, the ISW shoals over a realistic gentle bathymetric slope and complex background stratification, background current profiles sampled in the South China Sea. As the ISW shoals, the resulting convective and shear instabilities and ISW-driven turbulence are investigated.