Gas–liquid bubbly flow in two-dimensional bubble columns is studied by numerical simulation. An Eulerian–Eulerian two-fluid model is used to describe the time-dependent motion of the liquid driven by small, spherical gas bubbles injected at the bottom of the columns. The simulations are able to capture the large scale structures as observed experimentally in the laboratory. The numerical results, which include the characteristics of large scale structures such as wave length and frequency, mean velocities, turbulence intensities and turbulence shear stress, are analized and compared with the experimental data of Lin et al. (1996) [A.I.Ch.E. J., 42, 301–318] and Mudde et al. (1997) [A.I.Ch.E. J., 43, 913–926] and show good agreement.