In this paper, we have presented a comprehensive analysis of the development of flow pattern in a bubble column reactor by the introduction of a population balance equation combined with the three-dimensional two-fluid model (Reynolds stress model). The multiple size group (MUSIG) model has been used to account for the nonuniform bubble size distribution in a gas-liquid mixture. The coalescence and breakage effects of the gas bubbles are modeled according to the coalescence by the random collision driven by turbulence and wake entrainment while for bubble breakage by the impact of turbulent eddies. Local radial distributions of the gas hold-up, Sauter mean bubble diameter, axial liquid velocity, turbulent kinetic energy, turbulent energy dissipation rate, and Reynolds stresses for superficial gas velocity of 20 mm/s are compared against experimental data in a bubble column reactor. The development of flow pattern were examined at six axial locations HID = 0.2, 1.4, 2.6, 3.9, 5.0, and 6.2. Good quantitative agreement with the experimental data is obtained with three different models (i.e., k-ε, RSM with constant bubble size, and RSM with population balance). The model prediction shows better agreement with the experimental data with population balance than constant bubble diameter predictions.