  We use large-scale simulations to investigate the morphology of reionization
during the final, overlap phase. Our method uses an efficient three-dimensional
smoothing technique which takes into account the finite mean free path due to
absorption systems, lambda, by only smoothing over scales R_s<lambda. The large
dynamic range of our calculations is necessary to resolve the neutral patches
left at the end of reionization within a representative volume; we find that
simulation volumes exceeding several hundred Mpc on a side are necessary in
order to properly model reionization when the neutral fraction is ~0.01-0.3.
Our results indicate a strong dependence of percolation morphology on a large
and uncertain region of model parameter space. The single most important
parameter is the mean free path to absorption systems, which serve as opaque
barriers to ionizing radiation. If these absorption systems were as abundant as
some realistic estimates indicate, the spatial structure of the overlap phase
is considerably more complex than previously predicted. In view of the lack of
constraints on the mean free path at the highest redshifts, current theories
that do not include absorption by Lyman-limit systems, and in particular
three-dimensional simulations, may underestimate the abundance of neutral
clouds at the end of reionization. This affects predictions for the 21 cm
signal associated with reionization, interpretation of absorption features in
quasar spectra at z ~5-6, the connection between reionization and the local
universe, and constraints on the patchiness and duration of reionization from
temperature fluctuations measured in the cosmic microwave background arising
from the kinetic Sunyaev-Zel'dovich effect.
