The dynamics of a single or a group of slippery spheres settling under gravity in a Newtonian fluid is studied numerically. We focus particularly on the effect of particle surface slip on the sedimentation behavior. The flows containing moving slippery spheres are solved by a three-dimensional lattice Boltzmann model where a kinetic boundary condition is used to handle the slip phenomenon at the curved particle surface. The method is first validated by simulating the slip flow in a cylindrical tube, and the no-slip flows around one and two spheres settling in a container. The hydrodynamic behaviors of one, two and multiple slippery spheres settling under gravity are then investigated. The results for a single sphere show that the surface slip makes the sphere fall faster than a no slip particle and the wall correction factor decreases as the level of particle-surface slip is increased, indicating a drag reduction caused by the slip condition. For two settling spheres, when the no-slip particle is placed below the slip one, the two spheres will enter into the kissing phase earlier; on the contrary, deploying the no-slip particle above the slip one, the DKT process does not occur beyond a critical slip level and initial gap distance. If the two spheres are both slippery, the settling dynamics are similar to the no-slip case, but the time duration of the kissing phase decreases. As for the sedimentation of multiple spheres, it is found that the initial geometric arrangement has a significant impact on the sedimentation behavior. In general, slippery spheres in a cluster will experience larger fluctuations in the vertical velocity and position in the accelerated-falling stage, and smaller fluctuations in the decelerated-falling stage. (C) 2017 Elsevier B.V. All rights reserved.