The Bubble Population Balance Equation is implemented both in the two-fluid model and in the Algebraic Slip Mixture Model. Two and three-dimensional numerical simulations of two-phase (air–water) and pseudo-two-phase (air–Therminol-Glass beads) transient flows are performed for laboratory scale bubble columns of two different diameters operated over a range of superficial gas velocities (8 to 20cm/s) in the churn turbulent regime. The ability of the simulations to predict liquid recirculation and gas radial holdup profile, which drives the recirculation, is illustrated. The simulations also capture the effects of increased superficial gas velocity and of surface tension on bubble size distribution in agreement with observations. At all conditions a uni-modal bubble size distribution is predicted.