<p>Bubble column reactors are highly attractive for a wide range of industrial applications due to the combination of a simple structure, low operating costs, and excellent mass and heat transfer coefficients. Optimizing operation and developing new applications requires a thorough understanding of the fluid dynamics within these devices. However, this is a challenging task due to complex phenomena occurring simultaneously in a two-phase flow inside a bubble column, especially under heterogeneous flow regimes where bubble interactions and turbulence intensify. In this context, this study aims to contribute to the comprehension of the impact of closure models for drag and lift coefficients and critical Weber number on the flow dynamics inside a bubble column reactor, using computational fluid dynamics (CFD) simulations coupled with a population balance equation (PBE). The results have shown that the drag models (Ishii–Zuber and Tomiyama) perform similarly at low velocities but overestimate gas holdup at higher velocities. Compared to results from the Tomiyama lift model, a constant negative coefficient improved gas distribution by concentrating bubbles in the column’s center. The effects of increasing the critical Weber number are more evident when the Tomiyama’s model is applied for the lift coefficient. This is due to the resulting changes in bubble diameter and, consequently, in their distribution within the column. These changes lead to values that approximate those obtained using a constant value for the lift coefficient. This study underscores the critical role of closure models selection in achieving accurate multiphase simulations.</p>