  Reionization is inhomogeneous for two reasons: the clumpiness of the
intergalactic medium (IGM) and clustering of the discrete ionizing sources.
While numerical simulations can in principle take both into account, they are
at present limited by small box sizes. On the other hand, analytic models have
only examined the limiting cases of a clumpy IGM (with uniform ionizing
emissivity) and clustered sources (embedded in a uniform IGM). Here, we present
an analytic model for the evolving topology of reionization that includes both
factors. At first, recombinations can be ignored and ionized bubbles grow
primarily through major mergers. As a result, reionization resembles
"punctuated equilibrium," with a series of well-separated sharp jumps in the
ionizing background. These features are local effects and do not reflect
similar jumps in the global ionized fraction. We then combine our bubble model
with a simple description of recombinations in the IGM. We show that the
bubbles stop growing when recombinations balance ionizations. If the IGM
density structure is similar to that at moderate redshifts, this limits the
bubble radii to ~20 comoving Mpc; however, it may be much smaller if the IGM is
significantly clumpier at higher redshifts. Once a bubble reaches saturation,
that region of the universe has for all intents and purposes entered the
"post-overlap" stage, so the overlap epoch actually has a finite width. This
picture naturally explains the substantial large-scale variation in
Lyman-series opacity along the lines of sight to the known z>6 quasars. More
quasar spectra will shed light on the transition between the "bubble-dominated"
topology characteristic of reionization and the "web-dominated" topology
characteristic of the later universe. [Abridged]
