Epitaxial growth is a promising strategy to produce high-quality graphene samples. At the same time, this method has great flexibility for industrial scale-up. To optimize growth protocols, it is essential to understand the underlying growth mechanisms. This is, however, very challenging, as the growth process is complicated and involves many elementary steps. Experimentally, atomic-scale in situ characterization methods are generally not feasible at the high temperature of graphene growth. Therefore, kinetics is the main experimental information to study growth mechanisms. Theoretically, first-principles calculations routinely provide atomic structures and energetics but have a stringent limit on the accessible spatial and time scales. Such gap between experiment and theory can be bridged by atomistic simulations using first-principles atomic details as input and providing the overall growth kinetics, which can be directly compared with experiment, as output. Typically, system-specific approximations should be applied to make such simulations computationally feasible. By feeding kinetic Monte Carlo (kMC) simulations with first-principles parameters, we can directly simulate the graphene growth process and thus understand the growth mechanisms. Our simulations suggest that the carbon dimer is the dominant feeding species in the epitaxial growth of graphene on both Cu(111) and Cu(100) surfaces, which enables us to understand why the reaction is diffusion limited on Cu(111) but attachment limited on Cu(100). When hydrogen is explicitly considered in the simulation, the central role hydrogen plays in graphene growth is revealed, which solves the long-standing puzzle into why H2 should be fed in the chemical vapor deposition of graphene. The simulation results can be directly compared with the experimental kinetic data, if available. Our kMC simulations reproduce the experimentally observed quintic-like behavior of graphene growth on Ir(111). By checking the simulation results, we find that such nonlinearity is caused by lattice mismatch, and the induced growth front inhomogeneity can be universally used to predict growth behaviors in other heteroepitaxial systems. Notably, although experimental kinetics usually gives useful insight into atomic mechanisms, it can sometimes be misleading. Such pitfalls can be avoided via atomistic simulations, as demonstrated in our study of the graphene etching process. Growth protocols can be designed theoretically with computational kinetic and mechanistic information. By contrasting the different activation energies involved in an atom-exchange-based carbon penetration process for monolayer and bilayer graphene, we propose a three-step strategy to grow high-quality bilayer graphene. Based on first-principles parameters, a kinetic pathway toward the high-density, ordered N doping of epitaxial graphene on Cu(111) using a C5NCl5 precursor is also identified. These studies demonstrate that atomistic simulations can unambiguously produce or reproduce the kinetic information on graphene growth, which is pivotal to understanding the growth mechanism and designing better growth protocols. A similar strategy can be used in growth mechanism studies of other two-dimensional atomic crystals.