We explore the cause of the solar cycle variabilities using a novel 3D Babcock–Leighton dynamo model. In this model, based on the toroidal flux at the base of the convection zone, bipolar magnetic regions (BMRs) are produced with statistical properties obtained from observed distributions. We find that a little quenching in BMR tilt is sufficient to stabilize the dynamo growth. The randomness and nonlinearity in the BMR emergences make the poloidal field unequal and cause some variability in the solar cycle. However, when observed scatter of BMR tilts around Joy’s law with a standard deviation of 15°, is considered, our model produces a variation in the solar cycle, including north-south asymmetry comparable to the observations. The morphology of magnetic fields closely resembles observations, in particular the surface radial field possesses a more mixed polarity field. Observed scatter also produces grand minima. In 11,650 years of simulation, 17 grand minima are detected and 11% of its time the model remained in these grand minima. When we double the tilt scatter, the model produces correct statistics of grand minima. Importantly, the dynamo continues even during grand minima with only a few BMRs, without requiring any additional alpha effect. The reason for this is the downward magnetic pumping which suppresses the diffusion of the magnetic flux across the surface. The magnetic pumping also helps to achieve 11-year magnetic cycle using the observed BMR flux distribution, even at the high diffusivity.