In bubbly flow systems, bubbles may break or coalesce due to bubble?bubble or bubble?fluid interactions in presence of turbulence. In general, the population balance equations (PBEs) need to be solved to model the bubble size distribution (BSD). Two methods for solving the PBEs, i.e., the classes method and the quadrature method of moments (QMOM), were compared in this work. In the classes method, BSD was represented through a finite number of bubble classes, and coalescence rates and breakup rates were transformed into birth and death rates for each class. QMOM solved the equations of lower-order moments of BSD, instead of tracking the representative bubble size. A three-dimensional simulation of two-phase flow (air?water) was performed for a cylindrical gas?liquid bubble column reactor operated at high superficial gas velocities, and the computational fluid dynamics (CFD) were coupled with PBEs. An Euler-Euler two fluid model approach with an RNG k?? model of turbulence model was used. The predictions of the classes method (20 bins), QMOM with four moments and QMOM with six moments were compared with the experimental data in literature. The results showed that both the classes method and QMOM can reasonably predict the time-averaged gas volume fraction, liquid velocity profiles, mean diameter and bubble size distribution. Nevertheless, compared with the classes method, QMOM can save 2~3 times of computational resources. And for the QMOM four moments suffices were used to accurately describe the evolution of the gas phase. In fact, the results found by using the QMOM with four moments and the QMOM with six moments were very similar. In addition, the continuous bubble size distribution was reconstructed by using the low-order moments of QMOM. The predicted BSD was similar to that of the classes method, demonstrating that QMOM was a more efficient method than classes method when the PBM was coupled with CFD simulations for bubbly flow systems.