Coupling of the Discrete Element Method with Computational Fluid Dynamics (DEM–CFD) is a widely used approach for modeling particle–fluid interactions.Although DEM–CFD focuses on particle–fluid interaction, the particle–particle contact behavior is usually modeled using a simple Kelvin–Voigt contact model which may not represent realistic interactions of particles in high viscosity fluids. This paper presents an implementation of a new user-defined contact model that accounts for the effects of lubrication of fluid between two approaching particles while maintaining all other DEM–CFD particle–fluid interaction phenomena. Theoretical model that yields a non-linear restitution coefficient for submerged particle collisions, which was developed by Davis et al. (J Fluid Mech 163:479–497, 1986), is implemented in a DEM–CFD code. In this model, the behavior of particles at a contact depends on fluid properties, particle velocities and distance between particle surfaces. When two particles approach each other in a fluid, their kinetic energy decreases gradually because of a lubrication effect associated with the thin fluid layer between the particles. Particle post-collision behavior is governed by a simplified elastic contact law. With lubrication, it is possible that particles are not able to rebound if the approaching velocity is completely damped by lubrication, and in this case the particles agglomerate in the fluid. Tangential surface friction-slip forces are activated as in the case of dry particle contact. The lubrication model represents an advanced submerged particle collision approach that permits improved accuracy when modeling problems with high particle concentrations in a fluid. An application of the new model is shown in a simple sediment transport problem.