The lack of remote observations and the sparsity of in-situ data in the vast volume of the magnetosphere has made it difficult to apply data assimilation techniques, commonly used in meteorological and ionospheric research. Despite these limitations, in recent years the application of machine learning methods has enabled global reconstructions of plasma parameters from sparse datasets. One such success is the reconstruction of the global magnetic field and its dynamical evolution during storms and substorms by data-mining (DM) decades of space magnetometer data and fitting a flexible analytical description of the magnetic field. This approach, most recently implemented as SST19 (Stephens et al., 2019), has successfully reproduced the thinning and stretching of the magnetotail during the substorm growth phase, the formation of a tail X-line and subsequent rapid dipolarization during the substorm expansion phase, and the development of the ring current and plasma pressure distribution during storms. Meanwhile, magnetohydrodynamic (MHD) simulations of the magnetosphere have also greatly advanced in recent years. For instance, MHD codes, an example being GAMERA (Zhang et al., 2019), have revealed the importance of mesoscale structures, such as bursty bulk flows, during storms and substorms. However, the physical mechanisms responsible for storms and substorms are largely dictated by particle kinetic effects, which cannot be represented in the fluid approximation used in MHD. For this reason, much work has been done to couple MHD codes to kinetic models. Here, we present another approach, termed Gray-Box modeling, whereby the DM-based global configuration of the magnetic field from SST19 is used to guide GAMERA-MHD simulations of storms and substorms. We are pursuing this effort along three avenues. The first modifies the momentum equation to account for the formation of an embedded ion-scale current sheet, which “overstretches” the tail and is critical for the description of the substorm growth phase. The second initiates tail reconnection in the data-inferred location by including finite resistivity in those regions. The third, adjusts the entropy and flux-tube volume terms in GAMERA to account for the particle drift effects responsible for the formation of the ring current during storms.