In bubbly flows, the bubble size distribution dictates the interfacial area available for the interphase transfer processes and, therefore, understanding the behaviour and the average features of the bubble population is crucial for the prediction of these kinds of flows. In this work, by means of the STAR-CCM+ code, the S gamma population balance model is coupled with an Eulerian-Eulerian two-fluid approach and tested against data on upward bubbly pipe flows. The S gamma model, based on the moments of the bubble size distribution, tracks the evolution of the bubble sizes due to bubble break-up and bubble coalescence. Good accuracy for the average bubble diameter, the velocity and the void fraction radial profiles is achieved with a modified coalescence source. Numerical results show that better predictions are obtained when these flows are considered to be coalescence dominated, but, nevertheless, additional knowledge is required to progress in the development of coalescence and break-up models that include all the possible responsible mechanisms. In this regard, there is a requirement for experimental data that will allow validation of both the predicted bubble diameter distribution and the intensity of the turbulence in the continuous phase which has a significant impact on coalescence and break-up models. An advanced version of the model, that includes a Reynolds stress turbulence formulation and two groups of bubbles to account for the opposite behaviour of spherical bubbles, which accumulate close to the pipe wall, and cap bubbles, that migrate towards the pipe centre, is proposed. The Reynolds stress model is found to better handle the interactions between the turbulence and the interphase forces, and the use of only two bubble groups seems sufficient to describe the whole bubble spectrum and the bubbly flow regime up to the transition to slug flow. (C) 2016 Elsevier Ltd. All rights reserved.