Plunging liquid jets are commonly encountered in nature and are widely used in industrial applications (e.g., in waterfalls, waste-water treatment, the oxygenation of chemical liquids, etc.). Despite numerous experimental studies that have been devoted to this interesting problem, there have been very few two-phase flow simulations. The main difficulty is the lack of a quantitative subgrid model for the air entrainment process, which plays a critical role in this problem. In this paper, we present in detail a computational multiphase fluid dynamics (CMFD)-based approach for analyzing this problem. The main ingredients of this approach are a comprehensive subgrid air entrainment model that predicts both the rate and location of the air entrainment and a two-fluid transport model, in which bubbles of different sizes are modeled as a continuum fluid. Using this approach, a Reynolds-averaged Navier Stokes (RaNS) two-way coupled two-phase flow simulation of a plunging liquid jet with a diameter of 24 mm and a liquid jet velocity around 3.5 m/s was performed. We have analyzed the simulated void fraction and bubble count rate profiles at three different depths beneath the average free surface and compared them with experimental data in literature. We observed good agreement with data at all locations. In addition, some interesting phenomena on the different movements of bubbles with different sizes were observed and discussed.