Abstract: For nuclear reactor safety analysis, the one-dimensional two-fluid model equations are solved numerically with the first order upwind method because of its robust stability. In the present two-fluid model for horizontal stratified flow, surface tension is included because it makes the model well-posed. However, this is not done in industrial applications and numerical viscosity provides linear stabilization even when the model is ill-posed. It is now shown that numerical viscosity also provides nonlinear stabilization; meaning that the wave growth is bounded when the flow is unstable (e.g., in case of the Kelvin–Helmholtz instability). The formation of kinematic shocks in the presence of numerical viscosity provides the dissipation mechanism needed to stop the wave growth. However, numerical viscosity varies with the mesh size, which means that even though the unstable model is well-posed and meets von Neumann and nonlinear stabilization requirements, the solution does not converge for some short wavelengths greater than 2Δx. Furthermore, when the mesh size is large, significant artificial viscosity is added to the continuity and momentum equations and the model may become over-stabilized. A more scientific approach is proposed: to add a physical dissipation mechanism instead of numerical viscosity, i.e., Reynolds stresses. This approach has two advantages: first, the one-dimensional two-fluid model converges under severe dynamic conditions, such as the Kelvin–Helmholtz instability, and second, a higher order numerical method may be used instead of the first order upwind scheme without experiencing numerical excursions.