A simulation-based study of the variability of remnant polarization <inline-formula> <tex-math notation="LaTeX">({P}_{r}) </tex-math></inline-formula> in a multigranular 3-D ultrathin ferroelectric (FE) capacitor is presented in this article. The Poisson-Voronoi tessellation (PVT) algorithm is used for the nucleation of grains in the FE region, which corresponds to the physical growth mechanism. The PVT algorithm implemented in MATLAB is coupled with TCAD simulations, to trace the FE hysteresis loop. The impact of both, area and thickness scaling on the variability of <inline-formula> <tex-math notation="LaTeX">{P}_{r} </tex-math></inline-formula>, is considered. It is found that amount of variability in <inline-formula> <tex-math notation="LaTeX">{P}_{r} </tex-math></inline-formula> increases as FE thickness decreases. In addition, FE with the smaller surface area exhibits a higher variability in <inline-formula> <tex-math notation="LaTeX">{P}_{r} </tex-math></inline-formula> compared to a larger surface area FE capacitor. Furthermore, the impact of dielectric content in the FE grains is analyzed. It is seen that the dielectric grains cause a very large amount of variability in the FE hysteresis loop. An increase in the dielectric grains also leads to a loss in the retentivity of the hysteresis loop.