Modelling of gas-liquid bubbly flows is achieved by coupling a population balance equation with the three-dimensional, two-fluid, hydrodynamic model. For gas-liquid bubbly flows, an average bubble number density transport equation has been incorporated in the CFD code CFX 5.7 to describe the temporal and spatial evolution of the gas bubbles population. The coalescence and breakage effects of the gas bubbles are modeled. The coalescence by the random collision driven by turbulence and wake entrainment is considered, while for bubble breakage, the impact of turbulent eddies is considered. Local spatial variations of the gas volume fraction, interfacial area concentration, Sauter mean bubble diameter, and liquid velocity are compared against experimental data in a horizontal pipe, covering a range of gas (0.25 to 1.34m/s) and liquid (3.74 to 5.1m/s) superficial velocities and average volume fractions (4 to 21). The predicted local variations are in good agreement with the experimental measurements reported in the literature. Furthermore, the development of the flow pattern was examined at three different axial locations of L/D = 25, 148, and 253. The first location is close to the entrance region where the flow is still developing, while the second and the third represent nearly fully developed bubbly flow patterns.