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