The neutral mass density variation associated with the induced satellite drag force plays an important role in the Low Earth Orbit (LEO) spacecraft operations in the thermosphere, such as orbit maintenance, lifetime, and collision avoidance (Doornbos, 2012;Krauss et al., 2018;Zesta & Huang, 2016). Due to the impact of the solar-terrestrial energy input from the magnetosphere to the ionosphere-thermosphere system, the thermospheric density has intricate spatial and temporal variations. Consequently, predicting and simulating the neutral mass density changes and the related drag feedbacks to the spacecrafts are still a grand challenge for space weather research.Since the 1950s, several empirical thermospheric density models have been applied for orbit determination and prediction of LEO spacecraft, such as Jacchia model and Mass Spectrometer Incoherent Scatter (MSIS) model (Jacchia, 1971;Hedin, 1987Hedin, , 1991. However, the density discrepancies between the observations and these