A novel model is presented which uses mass and momentum equations based on population balances to describe the dispersed phase of bubble columns. A set of integro-differential non-linear equations constitutes the model, which can be solved directly using the least-squares methods. The implementation for these is presented. The model presents a degree of flexibility, as different density functions (number, mass and volume) and different internal coordinates (bubble diameter, volume and mass) can be used. A simplified 1D model of a bubble column is solved using the least-squares spectral element method. The results obtained for two different gas volume fluxes at the inlet are compared with experiments by several authors. The simulations show good agreement with the measured data. (C) 2011 Elsevier Ltd. All rights reserved.