Bubble column flows are encountered in a large number of industrial applications. In such flows gas sparged through the liquid rises in forms of bubbles of various sizes and provides the energy via interfacial momentum transfer for vigorous mixing of the liquid. Characteristic of bubble column flows are low (or zero) liquid superficial velocities, high superficial gas velocities and no moving mechanical parts. There is a need to predict the rate of liquid circulation and mixing driven by the unevenness in gas holdup cross sectional profiles and this is addressed by modeling of the two-phase fluid dynamics. This manuscript reviews the Computational Fluid Dynamics (CFD) modeling efforts and concepts in connection with such buoyancy driven gas-liquid flows in bubble columns. The published research during the past twenty years, with the main attention focused on the progress made during the past decade, is critically analyzed. The application capabilities of three different modeling approaches i.e., the Euler-Lagrange, the Euler-Euler and the Algebraic Slip Mixture Model, to computation of the dynamics of large-scale motion industrial installations have been reviewed. The issues associated with the interfacial momentum transfer and multiphase turbulence are discussed, and the need to incorporate more fundamental physics in such closures is stressed. It is realized that the interfacial and turbulence closures are still! challenges to surmount. Moreover, the improvements in the prediction capabilities of CFD models when using the bubble population balance coupled with the CFD computations are demonstrated. A need for further research to understand the bubble breakup and coalescence phenomena is emphasized.