When the exhaust plume from a descent engine impinges on the lunar surface, loose regolith can erode and become entrained into a high-velocity spray. These processes are simulated in this work by several integrated models: a hybrid continuum–kinetic solver for the gas flowfield, a coupled two-phase flow model for a polydisperse distribution of grain sizes, and a model for inelastic grain–grain collisions. The continuum regime is modeled with the data-parallel line relaxation code, and the kinetic modeling is done via the direct-simulation Monte Carlo method. Simulation results are first presented for a single-engine lander hovering at different altitudes. Surface stresses and the resulting dust erosion are compared to classical theory, and correction terms are introduced to improve agreement. The velocities of different-sized particles and particle mass fluxes are shown for different hovering altitudes. For a four-engine lander, there are multiple plume–plume and plume–surface interactions that result in complicated gas and dust flow structures. These interactions are discussed for three hovering altitudes, and the resulting dust sprays are compared to those of the equivalent thrust single-engine lander.