Abstract: The development of CFD-DEM is critical for investigating particle phenomena and their coupling with reactor transport. However, there continues to be considerable uncertainty in the selection of model parameters because of limitations in: (a) experimental measurements of multi-particle interactions, and (b) computational resources which have restricted most numerical studies to 2D simulations, in very small-scale systems (<50k particles) and/or to local sensitivity analysis. The focus of this study is to identify critical model parameters in 3D CFD-DEM simulations of fluidized beds through multivariate sensitivity analysis and quantify their impact on hydrodynamics. Towards this end, thirteen model parameters are considered and the sampling design matrix is constructed using the Morris-One-At-a-Time (MOAT) screening method. 3D CFD-DEM simulations with almost 170,000 glass bead particles (0.4 mm diameter) are conducted in a small rectangular pulsating fluidized bed, selected because of its repeatable bubbling patterns. Detailed bubble and particle dynamics data from 250+ simulations show that: (a) choosing exceedingly low normal spring stiffness has strong implications on particle velocities; (b) the impact of all contact dissipation parameters (normal restitution, friction and tangential damping) is tightly coupled and sensitivity to any one hinges on the choices of others; and (c) the stability of bubble patterns is contingent on their choices and almost-ideal as well as extremely dissipative systems exhibit no patterns. In addition, by investigating particle dynamics inside and around bubbles, we derive a working expression for the optimal choice of spring stiffness. Overall, this first-of-its-kind analysis provides important guidelines for CFD-DEM model parameter selection and the statistical framework developed here provides a robust strategy for the fundamental investigation of other particle-scale phenomena and simulation-based reactor design and optimization.