  We present results from a three-dimensional Babcock--Leighton dynamo model
that is sustained by the explicit emergence and dispersal of bipolar magnetic
regions (BMRs). On average, each BMR has a systematic tilt given by Joy's law.
Randomness and nonlinearity in the BMR emergence of our model produce variable
magnetic cycles. However, when we allow for a random scatter in the tilt angle
to mimic the observed departures from Joy's law, we find more variability in
the magnetic cycles. We find that the observed standard deviation in Joy's law
of $\sigma_\delta = 15^\circ$ produces a variability comparable to observed
solar cycle variability of $\sim $ 32%, as quantified by the sunspot number
maxima between 1755--2008. We also find that tilt angle scatter can promote
grand minima and grand maxima. The time spent in grand minima for
$\sigma_\delta = 15^\circ$ is somewhat less than that inferred for the Sun from
cosmogenic isotopes (about 9% compared to 17%). However, when we double the
tilt scatter to $\sigma_\delta = 30^\circ$, the simulation statistics are
comparable to the Sun ($\sim $18% of the time in grand minima and $\sim 10$% in
grand maxima). Though the Babcock--Leighton mechanism is the only source of
poloidal field, we find that our simulations always maintain magnetic cycles
even at large fluctuations in the tilt angle. We also demonstrate that tilt
quenching is a viable and efficient mechanism for dynamo saturation; a
suppression of the tilt by only 1-2$^\circ$ is sufficient to limit the dynamo
growth. Thus, any potential observational signatures of tilt quenching in the
Sun may be subtle.
