A 3D time dependent numerical study was performed to study the gas–liquid–solid three-phase flow dynamics in bubble columns by employing the Eulerian–Eulerian–Eulerian three-fluid approach. The three-phase bubble column of Gandhi et al. (Powder Technol. 1999, 103 (2), 80–94), Rampure et al. (Can. J. Chem. Eng. 2003, 81 (3–4), 692–706) and our previous study (Powder Technol. 2014, 260, 27–35) were studied. A mathematical model of a gas–liquid–solid three-phase flow was built. The effect of numerical schemes (wall boundary conditions, momentum discretization schemes and time steps) was discussed. CFD simulations were performed to study the sensitivity of the interphase drag models (different liquid–solid, gas–solid and gas–liquid drag models), and the predictions were compared with the experiments. The results showed that no-slip conditions for the wall boundary conditions, 0.001s for the calculation time step and second-order Upwind for momentum discretization, had better prediction results. The appropriate interphase drag forces to describe the three phase flow found in this study were the Zhang–Vanderheyden model, which was used as a gas–liquid drag model, Schiller–Naumann model, used as a liquid–solid drag model, and gas–solid drag force, which was not considered. The appropriate numerical schemes and interphase drag models were utilized to simulate the hydrodynamic parameters (time-averaged gas holdup, solid holdup, liquid axial velocity) in a gas–liquid–solid bubble column. The effects of superficial gas velocity, particle volume fraction, particle size and density were analysed and discussed. Four flow regimes of a gas–liquid–solid bubble column were simulated; the CFD results were compared with the experimental flow structures, and it was found that the bubble coalescence regime was better depicted. Superficial gas velocity had a large effect on the gas holdup of the bed, and the effect of solid volume fraction (Vs=0.03–0.30) and particle size (dp=75μm–270μm) on the distributions of time-averaged solid holdup and liquid axial velocity was greater than that of particle density (ρp=2500kg/m3–4800kg/m3). When particle size dp≥150μm and solid volume fraction Vs≥0.09, the hydrodynamic parameters had a strong dependency on dp and Vs. The larger the values of the Vs, dp and ρp were, the larger the axial solid concentration gradient was.