An adaptive numerical method is employed to simulate shear instabilities in open-ocean internal solitary-like gravity waves. The method is second-order accurate, employs block-structured adaptive mesh refinement, solves the incompressible Navier–Stokes equations and allows for the simulation of all of the length scales of interest by dynamically tracking important regions with recursively-nested finer grids. Two-dimensional simulations are performed over a range of parameters, which allows us to assess the conditions under which the shear instabilities in the waves occur, including a method to evaluate the critical Richardson number for instability based on the bulk wave parameters. The results show that although the minimum Richardson number is well below the canonical value of 1/4 in all simulations, this value is not a sufficient condition for instability, but instead a much lower Richardson number of 0.1 is required. When the Richardson number falls below 0.1, shear instabilities develop and grow into two-dimensional billows of the Kelvin–Helmholtz type. A linear stability analysis with the Taylor–Goldstein equation indicates that an alternate criterion for instability is given by σi Tw > 5, where σi is the growth rate of the instability averaged over Tw, the period in which parcels of fluid are subjected to a Richardson number of less than 1/4. A third criterion for instability requires that Lw/L > 0.86, where Lw is half the length of the region in which the Richardson number falls below 1/4 and L is the solitary wave half-width. An eigendecomposition of the rate-of-strain tensor is performed to show that the pycnocline thickness increases within the wave because of pycnocline-normal strain and not because of diffusion, which has important ramifications for stability. A three-dimensional simulation indicates that the primary instability is two-dimensional and that secondary, three-dimensional instabilities occur thereafter and lead to strong dissipation and mixing.