The final stages of cosmic reionization (EndEoR) are expected to be strongly regulated by the residual neutral hydrogen in the already ionized regions of the Universe. Its presence limits the mean distance that ionizing photons can travel and hence the extent of the regions that sources of ionizing photons can affect. The structures containing most of this residual neutral hydrogen are typically unresolved in large-scale simulations of reionization. Here, we investigate and compare a range of approaches for including the effect of these small-scale absorbers, also known as Lyman limit systems (LLSs), in such simulations. We evaluate the impact of these different approaches on the reionization history, the evolution of the ultraviolet background, and its fluctuations. We also compare to observational results on the distribution of Lyman-α opacity towards the EndEoR and the measured mean free path of ionizing photons. We further consider their effect on the 21-cm power spectrum. We find that although each of the different approaches can match some of the observed probes of the final stages of reionization, only the use of a redshift-dependent and position-dependent LLS model is able to reproduce all of them. We therefore recommend that large-scale reionization simulations, which aim to describe both the state of the ionized and neutral intergalactic medium, use such an approach, although the other, simpler approaches are applicable depending on the science goal of the simulation.