  By depleting the ionizing photon budget available to expand cosmic HII
regions, recombining systems (or Lyman limit systems) can have a large impact
during (and following) cosmic reionization. Unfortunately, directly resolving
such structures in large-scale reionization simulations is computationally
impractical. Instead, here we implement a sub-grid prescription for tracking
inhomogeneous recombinations in the intergalactic medium. Building on previous
work parameterizing photo-heating feedback on star-formation, we present
large-scale, semi-numeric reionization simulations which self-consistently
track the local (sub-grid) evolution of both sources and sinks of ionizing
photons. Our simple, single-parameter model naturally results in both an
extended reionization and a modest, slowly-evolving emissivity, consistent with
observations. Recombinations are instrumental in slowing the growth of large
HII regions, and damping the rapid rise of the ionizing background in the late
stages of (and following) reionization. As a result, typical HII regions are
smaller by factors of $\sim 2-3$ throughout reionization. The large-scale
($k\lesssim 0.2\text{ Mpc$^{-1}$}$) ionization power spectrum is suppressed by
factors of $\gtrsim 2-3$ in the second half of reionization. Therefore properly
modeling recombinations is important in interpreting virtually all reionization
observables, including upcoming interferometry with the redshifted 21 cm line.
Consistent with previous works, we find the clumping factor of ionized gas to
be $C_{\rm HII}\sim 4$ at the end of reionization.
