This thesis presents two important contributions to the modeling of entrainment of air bubbles in water, with focus on ship hydrodynamics applications. The first contribution consists of a general framework for modeling turbulent air entrainment. The framework attempts to describe the evolution of bubbles from their formation at the free surface, size distribution changes due to breakup and coalescence, and rise due to buoyancy. This proposed framework describes the complex entrainment process as a series of simpler mechanisms which can be modeled independently. For each mechanism a simple but mechanistic model is developed to provide closure while leaving the door open for future improvements. These unique characteristics enable the entrainment model to be used in general problems while still producing results at least as good as the few other available models. The massive entrainment of air that takes place around a ship leads to very high void fractions and accumulation of bubbles against the hull, particularly underneath the flat regions of the hull and in low pressure regions near appendages. These processes also pose challenges for two phase solvers. As a second contribution in this thesis, numerical algorithms for two phase flows are developed to eliminate the numerical instabilities normally occurring at high void fractions or large void fraction gradients. A hybrid method to improve pressure-velocity coupling for collocated grids is introduced, which keeps advantages typical of staggered grids in mass conservation and face flux computations. A new two phase coupling strategy is developed to guarantee stability at high void fraction. The balanced force method is extended to general curvilinear grids to suppress spurious velocities. The overall methodology provides strong coupling among pressure, velocity and void fraction, while avoiding numerical instability, and works for free-surface flows on dynamic overset grids. The proposed numerical schemes are tested for 1D and 2D cases. It is shown that the two phase solver is stable and efficient, even under extreme cases. Good mass conservation properties for multigroup simulations are also demonstrated. The air entrainment model is tested for a 2D wave breaking case and compared with extensive experimental data. The results show good predictions for entrainment location and two-phase properties. Full scale simulations for Athena R/V are performed using the same modeling constants obtained for the 2D wave breaking case. A grid study is also carried out to evaluate grid convergence properties of the model. While the model can predict well experimental data at full scale for the ship, it also shows dramatic improvements respect to previous entrainment models by converging in grid and not needing to re-evaluate the model constants for each new application. The high-speed Kann boat is also simulated at full scale, showing encouraging results for a preliminary entrainment model for aeration due to impact. The proposed numerical schemes are proved stable and robust in high Reynolds number flows with complex relevant geometries. In addition, these full scale simulations also identify modeling and numerical issues for future improvements.