High-amplitude internal solitary waves (ISWs) of depression shoaling over gentle slopes are examined using high-resolution numerical simulations. The computational tool consists of a three-dimensional nonlinear, non-hydrostatic hybrid spectral-element/Fourier method, online coupled with a high-accuracy Lagrangian particle tracker. The simulations are conducted over a long, realistic oceanic domain that follows a depth-dependent bathymetric transect of the South China Sea (SCS). The upslope propagation of the ISWs is aligned normal to the isobaths, while the fixed-in-time background stratification and current velocity are prescribed from in-situ measurements.At a critical depth determined by the wave amplitude, the ISW becomes convectively unstable (Umax > C), initiating baroclinically driven vorticity generation at the wave rear and the formation of a counter-rotating spanwise vortex pair. This structure defines the early primary two-dimensional instability. The three-dimensional evolution is triggered through a bypass transition to turbulence, and the flow structure becomes characterized by rising streamwise counter-rotating vortices and the formation of a turbulent subsurface recirculating core. As shoaling progresses and the ISW propagates over the steepest portion of the simulated transect, a vigorous plunge fills the wave interior with dense fluid originating from the rear shoulder. Nevertheless, the waveform remains remarkably symmetric due to the balance between nonlinear steepening and physical dispersion. At this stage, a gravity-current-like front slowly advances in a wave- following frame toward the leading edge of the ISW, expanding the recirculating core boundaries and producing enhanced turbulent mixing. Subsequent evolution indicates a significant reduction of the local stratification within the wave core, generating regions of sufficiently low Richardson number that trigger shear instabilities. These instabilities produce Kelvin–Helmholtz (KH) billows that disturb the trailing edge of the wave and drive a wake that modifies the background stratification. The present work gradually reveals key aspects of this dynamical evolution. Chapter 2 presents three-dimensional turbulence-resolving simulations of three separate ISWs with deep-water amplitudes that only vary by 5%. Convective and shear instabilities are examined both qualitatively and through metrics based on the kinetic energy associated with finer-scale motions. The turbulence formation within the subsurface core and the subsequent KH-driven wake are analyzed as a function of the initial ISW amplitude. Chapter 3 presents a fo- cused analysis using two-dimensional simulations with Lagrangian particle- tracking capabilities. For this base state of the ISW, the entrainment and de- trainment of neutrally buoyant particles within the materially coherent core, as well as their horizontal transport and residence times, are quantified. A complementary metric to characterize the core size is provided through Lagrangian coherent structures. Chapter 4 employs massively parallel three-dimensional simulations to investigate the sensitivity of the transition to turbulence to the domain width. Two additional turbulence-resolving simulations of the larger ISW examined in Chapter 2 are conducted in appreciably wider domains compared to the narrower cases of that chapter. An unrestricted development of streamwise vortices, producing lateral overturns not accommodated in the narrow domain due to width limitations, is observed. A quantification of spanwise length scales associated with energy-containing eddies is also presented. In addition, the larger number of such structures present leads to increased lateral variability and elevated kinetic energy inside and behind the wave, associated with velocity fluctuations in the wider domains. Finally, Chapter 5 summarizes the key findings and outlines directions for future investigations.