Bubble size distribution resulted from bubble coalescence and breakup has significant effects on the effective interphase drag in gas-solid bubbling fluidized beds, which was however not taken into account in previous coarse-grid simulations. In this study, the population balance model (PBM) was used to describe the dynamic evolution of gas bubbles in gas-solid bubbling fluidization, wherein the coalescence and breakup kernels were derived by considering the effects of bubble velocity difference, wake capture and bubble instability. An improved energy-minimization multi-scale (EMMS) bubbling model was developed to calculate the effective interphase drag in sub-grid scale. By incorporating both the PBM and the EMMS drag into an Eulerian continuum model, a CFD-PBM-EMMS coupled scheme was further proposed to predict the hydrodynamics of bubbling fluidized beds. The scheme was validated through comparison of simulated results with experimental data. The bed expansion characteristics and the lateral profiles of solids velocities were reasonably predicted at acceptable computational cost. Satisfactory agreement was also achieved between the measured and simulated bubble size distributions. The proposed sub-grid drag model and coupled simulation scheme can facilitate capturing the salient features of the hydrodynamics of gas-solid bubbling fluidized beds.