\u3cp\u3eWe present an efficient and highly scalable solver for direct numerical simulation (DNS) of dispersed gas-liquid flow, containing a large number of deformable bubbles. We apply this to O(10\u3csup\u3e4\u3c/sup\u3e) bubbles in a turbulent flow. This was accomplished by a well-considered combination of state-of-the-art numerical methods, as well as fast and scalable numerical algorithms that have their origin in single-phase and two-phase flows. The features key-elements of the algorithm of curvature computation, bubble collision and variable-coefficient Poisson equation solver are presented. Resolution requirements for an accurate advection of a rising bubble have been established by comparison with the (Hysing, 2009) benchmark. The Generalized Height Function (GHF) method, was adopted for the curvature computations. We observed agreement with theoretical convergence rates, as well as a significant reduction of spurious velocities when employing GHF, compared to a finite difference approach. The main interaction mechanisms between bubbles and walls of the domain were analyzed, showing second order convergence of the underlying numerical methods. A detailed analysis of the parallel performance of the Navier–Stokes (NS) solver and the solver for the gas volume fraction was carried out. The analysis revealed close to linear scaling up to ≈ 18000 cores on a computational grid with 1 billion cells for the NS solver, and an ideal scaling of the gas volume fraction solver up to O(10\u3csup\u3e3\u3c/sup\u3e) bubbles. Beyond that, acceptable overhead for up to O(10\u3csup\u3e4\u3c/sup\u3e) bubbles was found. A simulation of a downflow configuration of a turbulent channel loaded with a total of 10000 bubbles illustrates the computational capabilities. First and second order statistics of the velocity field were computed as well as the profiles of the average gas volume fraction field in the statistically steady state. In the case considered, a 47% increase of the wall shear stress was observed, brought about by turbulence modification arising from the embedded bubbles. The interaction between turbulence and bubbles at high volume fraction resulted in a strong attenuation of the rms of the velocity fluctuations in a wide region of the core of the channel.\u3c/p\u3