Modeling the late in-spiral and merger of supermassive black holes is central to understanding accretion processes and the conditions under which electromagnetic emission accompanies gravitational waves. In this regime, magneto-hydrodynamic numerical relativity simulations are needed to investigate how electromagnetic signatures correlate with black hole spins, mass ratios, and the gaseous environment. I will review the current status of simulations of binary black holes mergers in astrophysical environments.