Molecular dynamics simulations of plate-impact experiments were conducted to study shock propagation in boron carbide along two crystal orientations corresponding to the highest and lowest elastic moduli. While an elastic precursor wave was observed along both orientations upon shock loading, high intensity shocks also generated one and two trailing plastic waves in the former and latter orientations, respectively. These plastic waves were characterized by temperature rise, with accompanying loss of shear strength. The shock-velocity particle-velocity relationship for boron carbide, developed by measuring the distance travelled by these waves at different shock intensities, showed key similarities to experimental data. The axial stress-strain response through the loading and release stages were also determined, from which the shock Hugoniot curves were constructed. While the two orientations exhibited considerable anisotropy at low intensity shocks, high intensities revealed remarkable similarities in their Hugoniot states and release isentropes. This convergence in response was reasoned to arise from the loss of shear strength in boron carbide, which is in turn arose from the significant melting of the crystal due to conversion of plastic work into heat under high intensity shocks.