We report a parallel algorithm applicable to a Euler–Lagrange model embedding four-way coupling. The model describing the dispersed phase dynamics accounts for bubble–bubble collisions and is parallelized using a mirror domain technique while the pressure Poisson equation for the continuous phase is solved using a domain decomposition technique implemented in the PETSc library [S. Balay, K. Buschelman, W.D. Gropp, D. Kaushik, M.G. Knepley, L.C. McInnes, B.F. Smith, H. Zhang, PETSc Web page: http://www.mcs.anl.gov/petsc, 2001]. The parallel algorithm is verified and it is found that it gives the same results for both phases as compared to the serial algorithm. Furthermore the algorithm shows good scalability up to 32 processors. Using the proposed method, a homogeneous bubbly flow in a laboratory scale bubble column can be simulated at very high gas hold-up (37%) while consuming a reasonable amount of calculation wall time.