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 ntunerical 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.001 s 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 (V-s = 0.03-0.30) and particle size (d(p) = 75 mu m-270 mu m) on the distributions of time-averaged solid holdup and liquid axial velocity was greater than that of particle density (rho(p) = 2500 kg/m(3)-4800 kg/m(3)). When particle size d(p) >= 150 mu m and solid volume fraction V-s >= 0.09, the hydrodynamic parameters had a strong dependency on d(p) and V-s. The larger the values of the V-s, d(p) and pp were, the larger the axial solid concentration gradient was. (C) 2015 Elsevier B.V. All rights reserved.