Introduction and MotivationMagnetized, collisionless shocks are ubiquitous phenomena in the heliosphere and in astrophysical plasma. Oblique or quasi-perpendicular interplanetary shocks are associated with the most energetic Coronal Mass Ejections (CMEs) and have super-Alfvénic speeds (e.g., Bale et al., 2005;Krasnoselskikh et al., 2013), with high Mach number shocks potentially most geoeffective. The Wind spacecraft, launched in 1994, provides continuous solar wind measurements on a halo orbit around the L1 Lagrange point (the distance of about 0.01 AU from the Earth toward the Sun where the spacecraft's motion around the Sun matches that of the Earth), that include magnetic field vector, and plasma properties. The Deep Space Climate Observatory (DSCOVR), launched in 2015, also around L1, in addition to Earth observing instruments, carries the Plasma-Magnetometer instruments that provides high-cadence magnetic field and plasma in-situ data. Due to their strategic location near L1, the satellites provide early warning for solar activity that propagates toward the Earth, and ample detailed heliospheric shock data. However, these single-point observations cannot provide the complete picture of the multi-dimensional multi-ion shock structure and dynamics-a gap that can be possibly filled by numerical modeling.High Mach number M (where M indicates the fast magnetosonic speed Mach number) shocks in heliospheric plasma such as the solar wind (SW) and Earths' magnetosphere were modeled with 2.5D hybrid models in the past in order to investigate ion-kinetic properties of these shocks. In the hybrid model ions are modeled as large ensemble of particles embedded in the magnetic field and electrons are fluid. Previous models focused primarily on electron-proton plasma (e.g.,