Dense-phase pneumatic transportation of bulk materials in the form of slug flow has become a very important technology in industry. In order to understand the underlying mechanisms of slug flow, this paper presents a numerical study of slug flow in horizontal pneumatic conveying by means of discrete particle simulation. To be computationally efficient, the motion of discrete particles is three-dimensional and the flow of continuum gas is two-dimensional, and periodic boundary conditions are applied to both gas and solid phases horizontally. The proposed numerical model is qualitatively verified by comparing the calculated and measured results in terms of particle flow pattern and gas pressure drop. Then the influence of operational conditions such as gas and solid flowrates on slug behavior is numerically investigated. It is shown that slug velocity linearly increases with gas flowrate but is not sensitive to solid flowrate, and slug length increases with both gas and solid flowrates. The results qualitatively agree with the experimental observations. Finally, forces governing the gas-solid flow are analyzed on different length scales. It is shown that the movement of a slug is macroscopically controlled by the axial particle-fluid and particle-wall interactions, whereas the particleparticle interaction microscopically causes a slug to sweep up particles in a settled layer. The magnitudes of these interaction forces increase with gas and solid flowrates.