The behavior of air bubble clusters rising in water and the induced flow field are numerically studied using a three-dimensional two-way coupling algorithm based on a vortex-in-cell (VIC) method. In this method, vortex elements are convected in the Lagrangian frame and the liquid velocity field is solved from the Poisson equation of potential on the Eulerian grid. Two-way coupling is implemented by introducing a vorticity source term induced by the gradient of void fraction. Present simulation results are favorably compared with the measured results of bubble plume, which verifies the validity of the proposed VIC method. The rising of a single bubble cluster as well as two tandem bubble clusters are simulated. The mechanism of the aggregation effect in the rising process of bubble cluster is revealed and the transient processes of the generation, rising, strengthening, and separation of a vortex ring structure with bubble clusters are illustrated and analyzed in detail. Due to the aggregation, the average rising velocity increases with void fraction and is larger than the terminal rising velocity of single bubble. For the two tandem bubble cluster cases, the aggregation effect is stronger for smaller initial cluster distance, and both the strength of the induced vortex structure and the average bubble rising velocity are larger. For the 20 mm cluster distance case, the peak velocity of the lower cluster is about 2.7 times that of the terminal velocity of the single bubble and the peak average velocity of two clusters is about 2 times larger. While for the 30 mm cluster distance case, both the peak velocity of the lower cluster and two clusters are about 1.7 times that of the terminal velocity of the single bubble.