This study investigates two-phase air–water flow over a ten-step gabion stepped spillway with a 26.57° slope using Flow-3D, the k–ε turbulence model, Forchheimer equations for porous media, and an air entrainment solver. The numerical model was calibrated against Wüthrich and Chanson’s (2014) experimental data. A cell size of 0.012 m, comprising 1,477,056 cells, was selected as the optimal mesh, achieving an R² of 0.954 and an RMSE of 0.15. Nine simulations were performed using three discharge rates (0.018, 0.037, and 0.114 m³/s), representing nappe, transition, and skimming flow regimes, respectively, and three porosities (0.20, 0.375, and 0.50). Because porosity and effective aggregate diameter were varied simultaneously, permeability changes were attributed approximately 87% to porosity and 13% to aggregate diameter. Increasing permeability enhanced interflow and reduced overflow, weakening momentum flux, cavity vortices, and surface turbulent kinetic energy, which drives air entrainment. Increasing porosity from 0.375 to 0.50 reduced energy dissipation, mean velocity, mean air concentration, and turbulent kinetic energy by up to 21.8%, 22.5%, 14.6%, and 45.5%, respectively. These effects were strongest in nappe flow and weakest in skimming flow. Turbulent kinetic energy was the most sensitive parameter, while air concentration was the least sensitive. Finally, a new empirical correlation for unaerated flow length achieved a predictive R² of 0.994 and a mean absolute percentage error of 4.3%.