The production of graphene with open band gaps for the manufacturing of graphene-based electronic and optical devices requires synthesis methods to either control the number of layers to enrich AB-stacked bilayer or trilayer graphene or control the extent of functionalization of monolayer graphene. Solution-phase dispersion of graphene is promising for both methods to create printable electronics and nanocomposites. However, both methods face common challenges, including controlling the surface morphology, reducing the turbostratic layering, and enhancing the dispersion stability. To address these challenges at the molecular level, we successfully combined molecular simulations, theoretical modeling, and experimental measurements. First, we probed the surface structure and electrostatic potential of monolayer graphene dispersed in a sodium cholate (SC) surfactant aqueous solution, which exhibits 2D sheets partially covered with a monolayer of negatively charged cholate ions. Similar to the case of carbon nanotube functionalization, one may regulate the binding affinity of charged reactants for graphene functionalization by manipulating the surface morphology. Subsequently, we quantified the interactions between two graphene-surfactant assemblies by calculating the potential of mean force (PMF) between two surfactant-covered graphene sheets, which confirmed the existence of a metastable bilayer graphene structure due to the steric hindrance of the confined surfactant molecules. The traditional Derjaguin-Landau-Verwey-Overbeek (DLVO) theory was found to be adequate to explain the long-range electrostatic repulsions between the ionic surfactant-covered graphene sheets but was unable to account for the dominant, short-range steric hindrance imparted by the confined surfactant molecules. Interestingly, one faces a dilemma when using surfactants to disperse and stabilize graphene in aqueous solution: on the one hand, surfactants can stabilize graphene aqueous dispersions, but on the other hand, they prevent the formation of new AB-stacked bilayer and trilayer graphene resulting from the reaggregation process. Finally, the lifetime and time-dependent distribution of various graphene layer types were predicted using a kinetic model of colloid aggregation, and each graphene layer type was further decomposed into subtypes, including the AB-stacked species and various turbostratic species. The kinetic model of colloid aggregation developed here can serve as a useful tool to evaluate the quality of graphene dispersions for subsequent substrate-transferring or functionalization processes.