A coupled model of computational fluid dynamics and population balance model (CFD-PBM) scheme was developed to simulate bubble size distribution in gas-solid bubbling fluidized beds. The bubbling gas fluidized bed was divided into a discrete bubble phase and a dense emulsion including a pseudo solids continua modeled by the kinetic theory of granular flow (KTGF). A modified bubble coalescence and breakup model was proposed and implemented in the population balance model. A pseudo bubble-emulsion drag force model was derived based on the bubble size distribution and used for simulating the gas-solid momentum exchanges between gas bubbles and the solids belong to the emulsion phase. By taking into account the effects of bubble coalescence and breakup, the coupled CFD-PBM model was capable of predicting the hydrodynamic behavior of a bubbling gas-solid fluidized bed. A benchmark simulation showed good agreements between the computation results and the literature experimental data. Particularly, the predicted time-averaged bubble local hold-up maps, bubble size distributions and their evolution agreed well with the measurement data. Grid convergence simulation results demonstrated that the present model is able to predict the major hydrodynamic behavior of a 2D bubbling fluidized bed using coarse computational grids, which makes the present model a promising tool for applications in large-scale fluidized bed reactors. (C) 2019 Elsevier Ltd. All rights reserved.