In this study, we simulated slug flow using interFoam in OpenFOAM, which enables 3D two-phase flow analysis, and examined flow characteristics such as velocity and turbulent kinetic energy (TKE). Numerical simulations were conducted under a total of 12 superficial Reynolds number conditions. The simulations reproduced the phenomenon in which the elongated bubble region behind the liquid slug formed in a narrower and longer shape than that at the front, confirming that qualitative reproduction was feasible. As ReSL increased, slug length decreased and the number of slugs within the pipe increased. Under the ReSG = 1377 condition, the velocity ratio U/Um inside the liquid slug was below 1.2, whereas U/Um ≥ 1.5–2.0 appeared in the elongated bubble region. Under ReSG = 4320, the elongated bubble region exhibited U/Um ≥ 2.0, while the liquid slug showed a relatively low velocity ratio of approximately U/Um ≈ 1.0. TKE/Ud2 was highest primarily along the liquid–gas interface. Under the ReSG = 1377 condition, the liquid slug interior showed low TKE/Ud2 < 0.05, whereas values exceeding 2 occurred at the liquid–gas interface. For ReSG = 4320, TKE/Ud2 = 0.15–0.3 or higher appeared at the elongated bubble region, while the liquid slug interior showed TKE/Ud2 ≈ 0.01–0.1. For VL/V, under ReSG = 1377, the value increased with ReSL, and the temporal fluctuation between the maximum and minimum values was large. Under ReSG = 4320, VL/V also increased with ReSL, but the fluctuation was not pronounced. These results indicate that a three-dimensional numerical model for two-phase flow can reasonably reproduce the complex behavior and flow characteristics of slug flow. From these findings, future work involving the correlation analysis of numerical data and drone-based imaging data may contribute to identifying crack formation mechanisms in underground stormwater tunnels.