Gas-liquid bubbly flow in 2-D bubble columns was studied by numerical simulation. A Eulerian-Eulerian two-fluid model used describes the time-dependent motion of liquid driven by small, spherical gas bubbles injected at the bottom of the columns. Such equations, numerically implemented in this work were derived by Zhang and Prosperetti. A distinctive feature of this method is the derivation of the disperse-phase momentum equation by averaging the particle (here, the bubble) equation of motion directly, not the macroscopic equation for the particle phase. Both the time-averaged quantifies and dynamic characteristics of the macroscopic coherent structures agree with the experimental data of Lin et al. and Mudde et al. The comparison of simulated results with data demonstrates that this physical model and numerical approach can provide the key features of the time-dependent behavior of dispersed bubbly flows qualitatively with reasonable quantitative accuracy. Effects of the number of injectors, magnitude of bubble-induced viscosity, and various parameters in the interphase momentum exchange were also studied by simulating various cases and comparing with measurements. The applicability of different boundary conditions and the sensitivity to the mesh system used are also examined.