Reynolds-averaged Navier-Stokes (RANS) modeling is expected to deliver an ensemble-averaged result for the majority of turbulent flows. This could lead to the conclusion that multi-cycle internal combustion engine (ICE) simulations performed using RANS must exhibit a converging numerical solution after a certain number of consecutive cycles. However, for some engine configurations unsteady RANS simulations are not guaranteed to deliver an ensemble-averaged result.In this paper it is shown that, when using RANS modeling to simulate multiple engine cycles, the cycle-to-cycle variations (CCV) generated from different initial conditions at each cycle are not damped out even after a large number of cycles. A single-cylinder GDI research engine is simulated using RANS modeling and the numerical results for 20 consecutive engine cycles are evaluated for two specific operating conditions. One condition is characterized by stoichiometric operation and stable combustion (COVIMEP < 2%) and the other features dilute combustion resulting in misfire events and much higher COVIMEP. In both cases, multi-cycle RANS simulation results show cyclic variability.An in-depth analysis of the most significant physical and chemical quantities highlights that CCV are caused primarily by the variability of the in-cylinder flow. For stoichiometric combustion this does not greatly affect flame propagation, and the resulting fluctuations of the pressure traces are narrow. However, in the event of dilute combustion, the variability of the flow gains more importance, therefore delivering large cyclic fluctuations of in-cylinder pressure. Even though unsteady RANS simulations are not expected to predict as much variability in engine flows as LES, in particular for stable cases, they can still be used to capture typical combustion stability features in an internal combustion engine.This paper shows that the occurrence of CCV when using multi-cycle RANS modeling is not a numerical artifact or the effect of changes in the computational grid from cycle to cycle. Cold-flow analysis reveals that the variability of the simulated flow is intrinsic for the unsteady RANS approach and can only be damped out by increasing numerical viscosity, either by increasing the cell size or by adding up-winding. This is proved in a simplified non-ICE case such as the simulation of the vortex shredding from a cylinder in cross flow. Also, the effect of numerical viscosity on CCV is shown for multi-cycle engine simulations by using a simplified combustion model. Finally, it is shown that multi-cycle RANS can successfully predict the effect of changing specific engine parameters, such as the properties of the ignition system, on cyclic variability and combustion stability.