The shoaling of an internal solitary wave (ISW) of depression is explored in three-dimensions (3D) through high-accuracy, fully nonlinear, and nonhydrostatic simulations. Time-averaged background stratification and current profiles from field observations, along with measured bathymetry data from the South China Sea (SCS), are used. The computational approach is based on a high-resolution and high-accuracy deformed spectral multidomain penalty method incompressible flow solver. Recent field observations in the SCS indicate the presence of a convective instability followed by a subsurface recirculating core that may persist for more than tens of km and drive turbulent-induced mixing, estimated to be up to four orders of magnitude larger than that typically found in the ocean. The preceding convective instability occurs due to a sudden decrease in the wave propagation speed, below the maximum horizontal wave-induced velocity, and possible from the stretching of the near-surface vorticity layer of the baroclinic background current from the propagating ISW. Motivated by such observations, the present study examines the onset of the 3D convective instability that results in subsurface recirculating core formation, as the ISW propagates and shoals in the normal-to-isobath direction. A noise field is inserted in the wave-induced velocity and density field to force the evolution in 3D. The initial instability has a transitional structure that develops in the lateral direction. The evolution of the lateral instability and subsequent transition to turbulence in the breaking wave is compared with the wave structured observed in the field. As such, a preliminary understanding of the formation of recirculating cores in ISWs, the driver for subsequent turbulence, mixing, and particle transport in the interior is obtained.