The dynamics of cratering caused by a supersonic jet on a granular bed of solid particles is investigated using a high fidelity, two-phase model and numerical simulations. To model the gas phase, the Large Eddy Simulation (LES) approach is used and a kinetic theory based model is employed for the solid phase. Three-dimensional time-dependent simulations are conducted for the purpose of elucidating the morphological features of the crater and their evolution with time. Parametric variations of the initial conditions are performed to understand the effect of the initial solid volume fraction in the bed, of the coefficient of restitution, of the jet Mach number, of the particle diameter, and of the particle material density. Simulation snapshots of the craters are in qualitative agreement with images from the Mars Science Laboratory (MSL) mission and from Earth-based experiments. Analysis of the results shows that solid particle compaction at the crater base and side walls increases with increasing volume fraction of the undisturbed particle bed and reaches solid volume fractions close to 0.65. The coefficient of restitution is shown to have no effect on the large scale dynamics of the crater but affects the small scale features, particularly the shape of elongated bursts of solid particle clouds as well as radially aligned structures on the outer crater walls. A smaller jet Mach number results In a smaller crater characterized by ridge-like structures and walls at an angle of approximately 45 with respect to the undisturbed bed whereas larger Mach number jets create craters with walls approximately perpendicular to the undisturbed bed. A smaller particle size and a lighter particle material density also result in wider and deeper craters. The formation of all morphological structures is explained in detail and is related to the physical phenomena from which they originate. Power law curve fits are also presented for the crater outer diameter as well as for the depth. The mean and rms profiles of the solid volume fraction are computed. and it is found that the peak rms in the solid volume fraction does not occur at the surface but rather in the depth of the crater. In addition, the mean and rms of the particle momentum flux are also computed and discussed. Furthermore, the sources of particle momentum are evaluated, and the viscous drag as well as the gas pressure gradient are observed to be the prominent effects. The results show that the model developed is robust and is able to simulate craters due to a range of external conditions operating over particle beds having different characteristics. (C) 2017 Elsevier Ltd. All rights reserved.